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INTRODUCTION 

The  development  of  non-invasive,  biomedical  optical  imaging  from  frequency-domain  photon 
migration  (FDPM)  measurements  of  near-infrared  (NIR)  light  propagation  depends  upon  (i)  the 
measurements  of  optical  signals  on  the  boundary  of  tissues  and  (ii)  the  numerical  techniques  enabling  the 
reconstruction  of  interior  optical  properties  from  such  measurements.  From  the  mapping  of  interior  optical 
properties,  it  is  envisioned  that  diseased  tissues  can  be  identified  and  diagnosed  based  upon  the  differences  in 
absorption  and  scattering  properties.  Under  prior  NIH  support,  we  have  developed  instrumentation  for  rapid 
acquisition  of  FDPM  data.  Briefly,  FDPM  consists  of  launching  intensity-modulated  light  at  the  air-tissue 
interface  and  detecting  the  phase-delay  and  amplitude  attenuation  at  another  point  distant  from  the  incident 
point  source.  In  the  Purdue  Photon  Migration  Laboratory  (PPML),  we  have  developed  rapid  multi-pixel 
methods  for  acquiring  large  data  sets  of  phase-delay  and  amplitude  attenuation  across  a  tissue  surface  for  use 
in  an  inversion  algorithm  in  order  to  perform  image  reconstruction.  In  addition,  since  we  have  found  that  the 
endogenous  contrast  offered  by  absorption  and  scattering  may  be  insufficient  for  biomedical  imaging,  we 
have  invented  a  method  for  inducing  contrast  using  fluorescent  contrast  agents.  By  using  fluorescent 
contrast  agents,  we  have  shown  that  the  inverse  problem  may  be  better  posed  and  that  biodiagnostic 
information  can  be  obtained  from  assessing  the  fluorescent  decay  kinetics  within  the  tissue. 

BODY 


The  inverse  imaging  problem  is  crucial  to  the  development  of  biomedical  optical  imaging,  the  subject 
of  USAMRMC  support.  To  date  we  have  developed  an  algorithm  using  efficient  numerical  techniques  for 
the  image  reconstruction  for  imaging  tissues.  Our  algorithm  has  been  developed  in  two  stages,  consisting  of 
the  forward  and  inverse  imaging  problems.  The  forward  problem  involves  solution  of  the  optical  diffusion 
equations  for  excitation  and  emission  light  using  a  uniform  (and  coarse)  finite  element  and  finite  difference 
meshed.  We  have  refined  our  forward  solver  to  include  graded  meshing  and  three-dimensional  geometries. 
To  date  we  have  formulated  the  inverse  problem  in  a  number  of  ways: 

(1)  As  an  optimization  problem  in  which  parameter  updates  of  absorption  and  lifetime  are  separately 
determined  through  a  linear  Newton-Raphson  approach  (see  attached  manuscripts  by  Paithankar,  et  al., 
1998  and  Troy  and  Sevick,  1999); 

(2)  As  a  non-linear  optimization  problem  in  which  parameter  updates  of  absorption  and  lifetime  are 
separately  determined  through  a  Truncated  Newton  algorithm  using  efficient  finite  element  formulation 
(see  manuscripts  by  Roy  and  Sevick-Muraca,  1999a,  1999b); 

(3)  As  a  linear  optimization  problem  in  which  parameter  updates  of  absorption  and  lifetime  are  determined 
simultaneously  and  in  which  the  forward  problem  is  posed  as  Born  scattering  (see  Lee  and  Sevick- 
Muraca,  1999); 

(4)  As  a  non-linear  optimization  using  Bayesian  conditioning  and  parameter  zonation  structure  to  obtain 
parameter  estimates  in  three  dimensional  space  (in  conjunction  with  Univ.  of  Vermont  and  NSF,  see 
Eppstein,  et  al.,  1999);  and 

(5)  As  a  simply  bound  constrained  optimization  for  two  and  three-dimensional  reconstruction  (see 
manuscript  by  Roy  and  Sevick,  1999). 

Manuscripts  are  included  in  the  Appendix. 

The  accuracy  and  performance  of  these  algorithms  have  been  investigated  using  2-D  “synthetic”  data.  We 
currently  are  extending  our  formulation  to  3-D  experimental  data  using  the  multi-pixel  imaging  system 
developed  under  NIH  support. 
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The  USARMC  research  is  significant  to  the  biomedical  engineering  field  in  that  it  seeks  to  develop  the 
mathematical  basis  to  answer  questions  of  optical  imaging  performance  in  realistic  3-D  geometries.  Since 
prior  work  on  the  inverse  imaging  problem  has  been  largely  confined  to  two  dimensional  systems,  the 
opportunities  to  employ  its  solution  using  actual  three-dimensional  data  has  been  limited.  The  development 
of  a  three-dimensional  reconstruction  algorithm  in  the  PPML  is  significant  since  the  experimental  acquisition 
of  FDPM  data  of  model  systems  is  readily  available  to  test  the  algorithm  in  upcoming  research. 

KEY  RESEARCH  ACCOMPLISHMENTS: 

(1)  Developed  a  forward  predictor  code  for  2-D  and  3-D  frequency-domain  photon  migration  based  upon 
finite  difference  and  finite  element  numerical  techniques. 

(2)  Developed  an  inversion  code  for  2-D  and  3-D  reconstruction  of  fluorescent  optical  properties  using 
frequency-domain  photon  migration  measurements. 

(3)  Tested  these  algorithms  successfully  on  synthetically  generated  data  sets  and  have  explored  their 
success  in  the  presence  of  experimental  measurement  noise. 

The  research  has  enabled  us  to  move  forward  under  different  support  to  build  a  fast  optical  imaging  device 
and  to  image  mammary  disease  in  vivo  within  a  canine  model.  Under  NSF  and  NIH  support,  we  have  tested 
these  algorithms  using  experimental  data.  While  this  was  a  wildly  successful  exploratory  grant,  our 
proposal  for  clinical  translational  research  for  sentinel  lymph  node  staging  was  not  accepted  for  further 
funding  at  the  white  paper  stage. 
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CONCLUSIONS: 

The  overall  goal  of  this  research  is  to  apply  recent  advances  in  frequency-domain  photon  migration  (FDPM) 
technologies  within  an  in  vivo  imaging  study  to  develop  the  basis  for  non-invasive,  image-guided  therapy 
and  diagnosis.  Over  the  past  nine  years,  research  to  (i)  engineer  rapid  measurement  devices  for  clinical 
measurements;  (ii)  develop  inversion  algorithms  for  accurate  image  reconstruction  and  spectroscopy;  and 
(iii)  discover  the  use  of  fluorescent  contrasting  agents  for  enhancing  image  quality  have  resulted  in  an 
imaging  system  now  ready  for  testing  in  vivo.  In  our  laboratory,  we  have  obtained  in  vivo  FDPM  images  of 
phase  and  modulation  of  the  detected  intensity  modulated  wave  at  100  MHz  that  was  created  from  re-emitted 
ICG  fluorescence  from  deep  “leaky”  vessels  feeding  reactive  and  involve  inguinal  canine  lymph  nodes.  The 
results  suggest  the  potential  use  of  the  near-infrared  technology  for  diagnosis  of  lymph  node  involvement  in 
melanoma  and  breast  cancer  patients.  Other  opportunities  for  FDPM  spectroscopic  imaging  also  exist.  By 
employing  intensity  modulated  light  at  the  peak  oxy-  and  deoxy-  hemoglobin  and  isosbestic  absorption 
wavelengths  (820,  754,  and  800  nm)  the  opportunities  to  provide  low  resolution  oxygenation  mapping  during 
the  course  of  radiation  and  chemo-therapies  in  order  to  guide  staging  is  afforded.  In  addition,  sentinel  lymph 
node  mapping  could  also  be  improved  with  NDR.  techniques,  if  our  veterinary  studies  are  translatable  to  the 
clinic.  In  conjunction  with  proposed  Indiana  University/Purdue  University  Institutional  (IUPUI)  studies  of 


FINAL  REPORT  “Fluorescence  lifetime  imaging  for  breast  cancer  detection  and  diagnosis” 
Eva  M.  Sevick-Muraca 
Page  8 


sentinel  lymph  node  mapping,  which  involves  the  injection  of  blue  dye  or  a  radionucleotide  for  regional 
lymph  node  status,  we  also  propose  to  adapt  fluorescence  FDPM  for  non-invasive  evaluation  lymph  node 
status  at  IUPUI.  We  propose  to  use  approved  contrast  agents  (ICG)  in  patients  with  metastatic  breast  disease 
for  assessing  the  lymph  node  status  and  multi-focality  of  the  disease.  Clearly  the  junxtaposition  of  FDPM 
imaging,  sentinel  mapping,  and  axillary  lymph  node  biopsies  will  provide  comparative  ROC’s  for  assessing 
the  potential  impact  of  optical  imaging  in  the  management  of  breast  cancer  and  melanoma.  Finally,  the 
development  of  fluorescent  contrast  agents  with  lifetimes  sensitive  to  local  biochemical  environments  (p02, 
pH,  etc.)  may  also  provide  function/structural  imaging  capabilities  that  can  aid  in  the  efficacious  choice  of 
therapy  and  its  staging. 
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1.  Introduction 

Over  the  past  decade,  several  investigators  have  ex¬ 
plored  the  use  of  exogenous  fluorescent  dyes  as  con¬ 
trast  agents  to  differentiate  diseased  and  normal 
tissues  from  noninvasive  or  endoscopic  optical  mea- 
surements.  The  diagnosis  of  bum  depth  following 
Indocyanine  Green  dye  administration1  and  the  de¬ 
marcation  of  neoplastic  tissues  following  intravascu¬ 
lar  porphyrin  dye  administration2-4  are  possible 
because  of  their  leakage  from  vessels  corrupted  by 
msult  or  disease.  The  concomitant  increase  in  fluo- 
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rescence  from  the  dye  can  be  detected  at  the  tissue 
surface  and  permits  the  detection  of  disease.  With 
the  development  of  dyes  and  photodynamic  agents 
that  excite  and  reemit  in  the  near-infrared  red  wave¬ 
length  regime,  the  noninvasive  detection  of  diseased 
tissues  located  deep  within  tissues  may  also  be  pos¬ 
sible  because  red  excitation  and  reemission  light  can 
travel  significant  distances  to  and  from  the  tissue-air 
interface. s  However,  a  long-standing  problem  has 
been  the  low  uptake  or  leakage  into  neoplastic  tis- 
sues,  providing  insufficient  contrast  for  the  detection 
of  diseased  tissues.  Although  targeted  delivery  mav 
improve  uptake  ratios  and  contrast  imparted  by 
these  otherwise  therapeutic  agents,  these  approaches 
have  been  elusive  even  in  conventional  imaging  mo¬ 
dalities  such  as  magnetic  resonance  imaging  (MRI) 
and  x-ray  computed  tomography.  Investigators 
have  sought  to  alleviate  the  optical  contrast  uptake 
problem  by  employing  dyes  that  fluoresce  differently 
m  diseased  tissues  than  in  normal  tissues.  'The  use 
of  agents  whose  reemission  characteristics  vary  with 
tissue  pH6’7  and  p028  may  not  only  provide  detection 
of  diseased  tissues  by  the  nature  of  differing  fluores- 


cent  properties  but  may  also  contain  diagnostic  infor¬ 
mation.  The  difficulty  lies  in  measuring  the 
multiply  scattered  reemitted  light  and  reconstructing 
an  image  that  differentiates  tissues  on  the  basis  of 
fluorescent  properties,  such  as  fluorescent  yield  and 
lifetime. 

In  this  study  we  investigate  the  combination  of 
photon  migration  imaging  techniques  with  fluores¬ 
cence  spectroscopy  to  examine  the  feasibility  of  fluo¬ 
rescence  lifetime  imaging  (FLI).9  FLI  may  be 
considered  an  optical  analog  to  MRI.  Whereas  MRI 
depends  on  monitoring  the  spatial  variation  in  the 
relaxation  time  of  spin-spin  states  to  provide  high- 
resolution  imaging  in  tissues,  FLI  depends  on  moni¬ 
toring  the  spatial  variation  in  the  fluorophore  lifetime 
and  yield.  The  complication  with  FLI  in  tissues 
arises  from  the  excitation  and  fluorescent  photon 
times  of  flight,  which  are  similar  in  magnitude  to  the 
fluorescence  lifetime.  Herein  we  combine  inverse 
imaging  techniques  (as  developed  previously  by 
Pogue  et  al .10)  and  fluorescence  lifetime  spectroscopy 
to  demonstrate  the  feasibility  of  FLI  by  using  com¬ 
putational  studies.  We  separately  introduce  the 
concepts  of  photon  migration  imaging  and  fluores¬ 
cence  spectroscopy  in  Section  2  and  then  describe  the 
combination  of  techniques  in  our  computational  ap¬ 
proach  to  FLI.  Our  reconstructed  images  show  the 
ability  to  image  both  fluorescence  lifetime  and  fluo¬ 
rescence  yield  independent  of  any  a  priori  informa¬ 
tion  and  point  to  the  development  of  fluorescent  dyes 
whose  reemission  characteristics  are  environmen¬ 
tally  sensitive  and  provide  contrast  for  the  optical 
detection  of  diseased  tissues. 

2.  Background:  Photon  Migration  Imaging  and 
Fluorescence  Spectroscopy 

The  ability  to  reconstruct  an  internal  map  of  absorp¬ 
tion  optical  properties  from  continuous  wave11  and 
absorption  and  scattering  properties  from  frequency- 
domain10-1213  measurements  has  been  successfully 
demonstrated  by  using  multiply  scattered  light  from 
random  media.  Continuous  wave  measurements 
consist  of  noninvasively  monitoring  the  attenuation 
of  light  as  a  function  of  position  around  the  periphery 
of  a  heterogeneous  tissue  phantom  model  in  response 
to  a  constant  intensity  source,  whereas  frequency- 
domain  measurements  consist  of  monitoring  the 
phase  and  amplitude  modulation  of  multiply  scat¬ 
tered  light  at  various  peripheral  positions  in  response 
to  an  incident  modulated  light  source.  Both  ap¬ 
proaches  employ  a  perturbation  analysis  whereby  the 
absorption  and  scattering  properties  at  each  pixel 
position  within  the  tissue  phantom  are  individually 
and  independently  adjusted  until  the  measurements 
of  reemitted  light  at  the  periphery  of  the  tissue  phan¬ 
tom  match  those  predicted  by  a  model  for  diffusive 
light  propagation  in  random  media.  These  studies 
show  the  potential  for  image  reconstruction  of  hidden 
tissue  heterogeneities  when  the  contrast  caused  by 
either  absorption  or  scattering  coefficients  is  as  low  as 
2:1. 13  However,  Troy  etal.14  conducted  in  vitro  mea¬ 
surements  of  normal  and  diseased  breast  tissues  and 


showed  that  the  endogenous  optical  contrast  between 
120  normal  and  diseased  breast  tissues  is  not  signif¬ 
icantly  different  for  consistent  detection  with  optical 
techniques.  These  measurements  were  conducted 
ex  vivo  and  consequently  may  underestimate  the  con¬ 
trast  available  in  vivo  caused  by  increased  tumor 
vasculature  and  absorption  caused  by  hemoglobin. 
Nonetheless,  these  results  suggest  the  need  for  opti¬ 
cal  contrast  agents. 

Previously,  we  experimentally  demonstrated  the 
increased  sensitivity  for  detecting  heterogeneities  on 
the  basis  of  fluorescence  as  opposed  to  absorption 
when  time-dependent  measurements  are  employed.15 
Indeed,  the  augmentation  of  optical  contrast  that  is 
due  to  the  lifetime  of  a  fluorophore  or  phosphorescent 
probe  has  been  recognized  by  Wu  et  al.16  in  tissue 
phantom  studies.  The  localization  of  a  fluorescing 
body  has  been  performed  by  various  researchers.17-19 
Upon  activation  into  a  higher  electronic  state  by  the 
absorption  of  light,  an  activated  fluorophore  may  un¬ 
dergo  nonradiative  decay  or  radiative  decay;  the  lat¬ 
ter  results  in  the  reemission  of  a  fluorescent  photon. 
The  yield  is  defined  by  the  fractional  number  of  flu¬ 
orescent  photons  reemitted  for  each  excitation  photon 
absorbed  or  the  fraction  of  decay  events  that  results 
in  the  emission  of  a  fluorescent  photon.  The  lifetime 
is  defined  as  the  mean  survival  time  of  the  activated 
fluorophore  or  the  mean  time  between  the  absorption 
of  an  excitation  photon  and  the  reemission  of  a  fluo¬ 
rescent  photon.  Because  the  stability  of  the  acti¬ 
vated  fluorophore  depends  on  its  environment,  both 
yield  and  lifetime  are  also  dependent  on  the  (bio) 
chemical  environment  in  which  the  fluorophore  re¬ 
sides.  Consequently,  the  lifetime  can  provide  con¬ 
trast  based  on  the  differences  in  biochemical 
environments  of  normal  and  diseased  tissues,  similar 
to  the  contrast  provided  by  relaxation  times  in  nu¬ 
clear  magnetic  resonance  imaging. 

The  use  of  lifetime  for  contrast  in  optical  imaging  in 
tissues  is  not  new.  In  studies  demonstrating  the 
noninvasive  differentiation  of  hematoporphyr in¬ 
laden  tumors  from  normal  tissue,  Cubeddu  et  al.20 
employed  the  comparatively  long  lifetime  (>15  ns)  of 
reemitted  fluorescence  caused  by  increased  porphy¬ 
rin  uptake  in  tumors  over  the  lifetime  of  endogenous 
compounds  (>7  ns)  in  normal  tissues  in  order  to  pro¬ 
vide  discrimination  of  the  two  tissue  types.  These 
microscopy  studies  involved  unscattered  light  that 
provided  correct  reemission  kinetics  of  measured  flu¬ 
orescence.  In  this  computational  study,  we  employ 
contrast  that  is  due  to  increased  uptake  as  well  as 
lifetime  in  order  to  differentiate  diseased  tissue  by 
using  multiply  scattered  light  reemitted  from  tissues 
that  involve  deeply  seated  tissue  volumes.  The  re¬ 
emission  kinetics  of  phosphorescent  probes  do  not 
permit  an  interrogation  deep  into  tissue  with  time- 
domain  approaches.21  Our  preliminary  frequency- 
domain  computations  (unpublished)  suggest  that  the 
simultaneous  determination  of  location  and  lifetime 
is  difficult  when  long-lived  phosphorescent  probes  are 
used.  In  Section  3  we  underscore  the  magnitude  of 
the  problem  by  describing  the  forward  imaging  prob- 
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lem  (i.e.,  prediction  of  the  frequency-domain  mea¬ 
surements  given  the  location,  optical  properties,  and 
fluorescent  properties  of  the  random  media)  as’ well 
as  our  adapted  approach  for  solving  the  inverse  prob¬ 
lem  (i.e.,  prediction  of  the  location,  optical  properties, 
and  fluorescent  properties  of  the  medium  from  mea¬ 
surements  of  frequency-domain  light  propagation). 

3.  Forward  Problem  and  the  Solution  Methodology 

The  spatial  and  temporal  transport  of  light  in  tissues 
or  multiply  scattering  media  can  be  accurately  de¬ 
scribed  by  the  diffusion  approximation  to  the  radiative 
transport  equation.  A  coupled  frequency-domain 
diffusion  equation  can  be  used  to  predict  the  excita¬ 
tion  and  emission  fluence  rates,  <1  \(r,  w)  and  <!>„,(/-,  w), 
at  any  location  r  within  a  tissue  phantom  by  Eqs.  (1) 
and  (2),  respectively.22-'24 

V  •  [£>,(r)Vc[>r(r,  w)]  -  [g„  (r)  +  iaj/cJ 

X  ‘lyr ,  w)  +  Sx(r,  w)  =  0,  (1 ) 

V  •  [Z),»VT,„(r,  w)]  -  [p.„  (r)  +  iw/cj 

X  <!>„,( r ,  to)  +  Sjr,  w)  =  0.  (2) 

The  source  term  for  excitation  light  St(r,  w)  is  due  to 
the  sinusoidally  modulated  light  at  frequency  w  = 
W,  where  fis  usually  in  the  megahertz  range.  The 
first  term  in  both  equations  represents  the  diffusive 
or  random-walk  transport  of  light  where  Dx  m  is  the 
optical  diffusion  coefficient,  i.e., 

=  [3(|x„rm  +  p^')]"1,  (3) 

and  and  gs'  are  the  absorption  and  isotropic  scat¬ 
tering  coefficients,  respectively.  The  optical  proper¬ 
ties  are  dependent  on  the  wavelength  of  light  and 
thus  are  different  for  the  excitation  (subscript  x)  and 
fluorescent  (subscript  m)  light.  The  total  absorption 
coefficient  at  the  excitation  wavelength,  ga  ,  is  due  to 
contributions  from  nonfluorescing  chromophores  as 
well  as  from  fluorescent  dye.  The  total  absorption 
coefficient  is  given  by  the  sum  of  absorption  coeffi¬ 
cients  that  are  due  to  nonfluorescing  chromophores, 
p-a  _,  and  fluorophores,  ga  ^ .  We  assume  that  the 
absorption  experienced  at  the  fluorescent  wavelength 
primarily  is  due  to  nonfluorescing  chromophores. 
The  velocity  of  light  in  tissue  is  cn  =  c/n,  where  n  is 
the  average  index  of  refraction.  The  source  term  for 
the  fluorescent  light  is  dependent  on  the  excitation 
light  fluence,  fl>x(r,  w),  and  is  given  by 

S„,(r,  to)  =  Hflu,_,„(r)<l>x(r ,  w)  — —  j  ^ .  (4) 

This  term  arises  from  the  Fourier  transform  of  the 
single-exponential  fluorescence  decay  term  (1/ 
T)exp(— 1/-[)  in  the  time  domain  following  an  incident 
pulse  of  excitation  light  where  r  is  the  fluorophore 
lifetime.  Here  t|  is  the  quantum  yield  and  the  ab¬ 
sorption  coefficient,  p,  ,  is  the  product  of  the  extinc¬ 
tion  coefficient,  loge  10,  and  the  concentration  of  the 
fluorophore  in  the  ground  state.  For  the  purposes  of 
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this  study,  the  combined  product,  tip.,,  ,  is  termed 
the  fluorescent  yield  and  is  proportional" to  the  gen¬ 
erated  fluorescence  fluence.  Note  that  multiexpo¬ 
nential  time  decay  can  also  be  handled  with  this 
procedure  by  a  simple  extension. 

In  the  source  term  for  fluorescent  light  given  in  Eq. 
(4),  p.„  changes  as  the  relative  fractions  of  the  flu¬ 
orophore  in  the  ground  and  excited  states  change. 
The  saturation  effects  would  have  to  be  handled  by 
taking  into  account  the  relative  amounts  of  fluoro¬ 
phore  in  the  ground  and  excited  states.  We  neglect 
the  saturation  effects  and  assume  single-exponential 
decay  kinetics  in  this  initial  study. 

Furthermore,  g  in  Eq.  (1)  is  the  sum  of  |in  and 
what  follows,  images  of  lifetime,  t,  and 
images  of  x|p,a^  are  obtained.  There  is  a  difficulty 
in  the  estimation  of  pn  that  is  given  by  p.„  =  p.„  + 

f-a,  because  the  explicit  values  of  p ,n  are  not 
known  (the  values  of  the  product  of  p  and  p.„  are 
known).  Because  the  major  contribution  to' (Ia  is 
from  and  not  p.„  we  have  used  an  approxi¬ 
mate  expression  for  p.0  ,  p.„t  =  p.a  _  +  pp.o  and  have 
chosen  t]  as  1  only  for  the  purposes  of  calculation  of 
(V,-  Our  approach  can  be  extended  to  permit  imag- 
lng  of  q,  p -O'  w,  p-ar  ,  and  t  from  measurements  con¬ 
ducted  at  both  excitation  and  fluorescent 
wavelengths.  Although  reconstruction  from  fluores¬ 
cent  measurements  provides  T]p.a  ,  excitation  mea¬ 
surements  can  provide  p.n  _  (in  the  absence  of  a 
fluorophore)  and  p.  +  ^  _  (in  the  presence  of  a 
fluorophore).  From  these  images,  maps  of  p  and 
can  be  obtained  in  principle. 

Both  Eqs.  (1)  and  (2)  are  linear  complex  elliptic 
equations  that  can  be  solved  as  boundary  value 
problems  for  complex  quantities  <t»c(r,  w)  and  <b„  (r, 
w).  We  employ  the  method  of  finite  differences  in 
which  we  place  a  grid  over  the  space  domain  and 
obtain  an  approximation  to  the  solution  at  each  grid 
Point,  j.  One  of  the  fastest  methods  to  solve  these 
linear  elliptic  boundary  value  problems  is  the  mul¬ 
tigrid  solution  (see  the  review  by  Fulton  et  al.25). 

In  the  procedure  for  the  multigrid  solution,  an  ini¬ 
tial  solution  is  obtained  quickly  for  coarse  grids  that 
is  then  further  refined  for  a  better  solution  for  finer 
grids.  This  is  an  involved  procedure  in  which  we 
have  elected  to  use  mudpack  routines.26  mudpack 
routines  are  flexible  and  allow  placement  of  sources 
either  at  the  surface  or  inside  the  phantom.  For 
the  equations  to  be  solved,  it  is  assumed  that  <l>m  x(r, 

(o)  —  0  on  the  tissue  surface,  which  is  known  as  the 
zero-fluence  boundary  condition.  This  is  imple¬ 
mented  by  assigning  the  absorption  coefficient  for 
both  excitation  and  fluorescent  light  at  all  the  grid 
points  in  the  square  grid  lying  outside  the  circular 
tissue  phantom  to  a  large  value.  The  source  was 
simulated  by  setting  the  value  of  <I>X  to  an  arbitrary 
complex  number  at  a  grid  point  on  the  surface 
where  the  source  is  located.  The  solution  of  Eqs. 

(1)  and  (2)  yields  a  complex  number  for  d>m  at  each 
grid  point,  j.  The  detected  signal  at  the  surface  is 
proportional  to  the  normal  component  of  the  gradi- 


Source  C 


Fig.  1.  Schematic  of  the  circular  simulated  tissue  phantom  inter¬ 
rogated  by  four  sources,  one  source  at  a  time.  Twenty  detectors 
are  located  on  the  periphery  with  uniform  separation.  A  circular 
object  that  mimics  a  hidden  diseased  tissue  volume,  located  within 
the  tissue  phantom,  is  also  shown. 


ent.  In  this  study,  for  evaluation  of  the  signal  at 
detector  i  located  on  the  surface,  we  used  the  value 
of  <Pm  at  an  internal  grid  point  closest  to  the  detec¬ 
tor.  This  is  reasonable  in  our  initial  approach  be¬ 
cause  the  normal  component  of  the  gradient  is 
proportional  to  just  inside  the  surface.10  The 
phase  lag,  dM,  and  the  log  of  the  ac  amplitude,  Mm, 
at  the  detectors  were  calculated  with  respect  to  the 
phase  and  the  ac  amplitude  of  the  source.  The 
computational  time  required  to  solve  for  the  fluo¬ 
rescent  fluence  on  a  SunSparc  10  for  a  17  X  17  grid 
is  0.04  s,  and  that  for  a  33  X  33  grid  is  0.165  s. 
Different  grid  sizes  were  used  in  the  testing  phase  of 
the  solution  of  the  equations,  and  the  solutions  were 
in  agreement  within  the  numerical  error  caused  by 
the  finite  grid  size. 

Before  attempting  to  solve  the  inverse  problem  for 
the  simulated  phantom  illustrated  in  Fig.  1,  we  must 
first  understand  the  effects  of  changing  the  fluores¬ 
cent  optical  properties  of  the  tissue  on  0m  and  Mm 
measured  at  a  detector  or  series  of  detectors  located 
on  its  surface.  Solutions  to  Eqs.  (1)  and  (2)  were 
obtained  in  two  dimensions  for  a  65  X  65  grid  cover¬ 
ing  a  100-mm-diameter  circular  tissue  phantom  with 
a  circular  embedded  heterogeneity  of  30  mm  diame¬ 
ter  and  located  at  the  center  of  the  tissue  phantom. 
The  simulated  measurements  of  fluorescent  phase 
shift  and  ac  amplitude  are  reported  for  20,  equally 
spaced,  circumferentially  located  detectors.  The 
modulation  frequency,  f,  was  set  equal  to  150  MHz. 
The  optical  properties  of  the  heterogeneity  and  the 
background  are  shown  in  Table  1. 


Fig.  2.  Log  plot  of  fluorescent  light  ac  amplitude  at  150  MHz  at 
various  detectors  for  a  heterogeneous  tissue  phantom  with  a  30- 
mm-diameter  object  located  at  the  center  of  the  phantom  for  var¬ 
ious  object  values.  Absorption  coefficient  r\\ia^  ^  is  1  X 

10 ~5  mm-1  for  the  background;  lifetimes  t  for  both  the  object  and 
the  background  are  1  ns. 


A.  Effect  of  Tun. a  ^  on  0m  and  Mm 

In  order  to  evaluate  the  influence  of  rj|xa  0m  and 
Mm  were  computed  at  each  detector  as  tKe  value  of 
Tl|jLa  m  in  the  heterogeneity  increased  from  1CT4 
mm-1  to  10"1  mm"1  and  as  T|pta  ^  in  the  background 
was  maintained  constant.  The  modulation  fre¬ 
quency  was  chosen  as  150  MHz.  The  lifetime,  t,  was 
set  equal  to  1  ns  for  both  the  object  and  the  back¬ 
ground  causing  contrast  resulting  solely  from  differ¬ 
ences  in  *n|jLa  n.  The  plots  of  0m  and  Mm  are  shown 
in  Figs.  2  and  3,  respectively,  for  one  source  located  at 
position  A.  The  curves  for  Mm  are  not  smooth  be¬ 
cause  the  circular  surface  is  approximated  by  joining 
discrete  grid  points  closest  to  the  circumference  by 
straight  lines.  Furthermore,  the  curves  for  both  0m 
and  Mm  are  symmetric  for  detectors  2  and  20,  3  and 
19,  4  and  18,  and  so  on  because  these  detector  pairs 
are  symmetric  with  respect  to  the  source.  The  de¬ 
tectors  near  the  source  are  unaffected  by  the  presence 
of  an  object  because  the  photon  sampling  volume  does 
not  include  the  region  occupied  by  the  heterogeneity. 
At  the  other  detector  locations,  the  fluorescent  ac 
amplitude  increases  as  m  of  the  heterogeneity 
increases.  As  of  the 'Heterogeneity  increases 

to  high  values,  the  ac  amplitude  approaches  an  upper 
limit  that  is  due  to  high  absorption  coefficient  jxa  , 
which  is  a  sum  of  \ia  _ ^  and  p,a  .  This  is  similar  to 
the  inner-filter  effect  in  which  the  high  absorption  in 
the  heterogeneity  shields  the  interior  of  the  hetero¬ 
geneity  from  excitation  photons.27  Figure  3  illus¬ 
trates  changes  in  fluorescent  phase  shift,  0,„,  as  a 


Table  1.  Optical  Properties  and  Experimental  Parameters  for  the  Forward  Problem 


M-nt 

’  Mr*,,, 

P-,/  or  mm' 

M'a,-. 

T)\LUr  ^  Background 
(mm-1) 

t  Background 

Frequency 

(mm1) 

( mm " 1 ) 

(mm"1) 

(mm”1) 

(ns) 

(MHz) 

P»,-.  +  Pa, . 

0.0 

1.0 

0.0 

1.0  x  10'5 

1.0 

150.0 
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Fig.  3.  Plot  of  fluorescent  light  phase  shia  at  150  MHz  at  various 
detectors  for  a  heterogeneous  tissue  phantom  with  a  30-mm- 
diameter  object  located  at  the  center  of  the  phantom  for  various 
object  T||x„t  values.  Absorption  coefficients  and  lifetimes  are  the 
same  as  in  Fig.  2. 
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Fig.  4.  Log  plot  of  fluorescent  light  ac  amplitude  at  150  MHz  at 
various  detectors  for  a  heterogeneous  tissue  phantom  with  a  30- 
mm-diameter  object  located  at  the  center  of  the  phantom  for  var¬ 
ious  object  t  values.  Background  lifetime  t  is  1  ns  and  absorption 
coefficients  are  1  X  10  r>  mm'1  for  the  background  and  1  X 

10  •*  mm  1  for  the  object. 


i=0.1  ns  — 
i=1  ns 
i=10ns  q 
t=100ns  - 
t=1000  ns 


i  a  , 


function  of  the  fluorescent  yield  that  is  due  to  the 
fluorophore,  r||xa  ,m-  As  T| jxa ^  is  increased  from 

10  nun  to  10  1  mm  1  (when  the  background 
value  is  10“'  mm-1),  little  change  in  phase  shift  oc¬ 
curs.  Because  the  phase  shift  of  reemitted  fluores¬ 
cence  in  dilute,  nonscattering  solutions  is  a  function 
of  lifetime  and  is  not  dependent  on  concentration  [0  — 
tan  (wt)],  one  might  expect  similar  trends  in  scat¬ 
tering  media.  Indeed,  as  shown  in  Fig.  3,  phase- 
shift  changes  caused  by  fluorophore  concentration 
differences  are  small.  Consequently,  the  changes  in 
phase  shift  with  increasing  ti fiQ  of  the  heterogene¬ 

ity  can  be  attributed  to  the  alteration  in  photon  mi¬ 
gration  as  a  result  of  the  presence  of  heterogeneities 
with  a  high  absorption  of  excitation  light.28  The 
presence  of  heterogeneities  with  high  absorption  re¬ 
duces  the  effective  path  length  of  photon  migration 
and  reduces  the  phase  shift,  though  this  reduction  is 
small  in  magnitude.  In  summary,  Mm  is  directly 
and  strongly  dependent  on  changes  in  -q|xa  of  a 
simulated  tissue  heterogeneity,  whereas  0m  is  indi¬ 
rectly  and  weakly  dependent  on  because  of 

changes  in  photon  migration. 

B.  Effect  of  r  on  Mm  and  0m 

In  order  to  evaluate  the  influence  of  t,  Mm  and  0 
were  calculated  at  each  detector  as  the  values  of  t  in 
the  heterogeneity  varied  from  10“ 1  ns  to  103  ns  and 
the  value  of  t  in  the  background  was  held  at  1  ns. 
The  modulation  frequency  was  at  150  MHz.  Back¬ 
ground  TUA„t  was  set  to  10“5  mm-1,  and  r)|xa  for 
the  object  was  set  to  10“3  mm-1.  As  shown  in  Fig.  4, 
the  detected  ac  amplitude  increases  as  r  decreases.' 
Because  Mm  of  reemitted  fluorescence  in  a  dilute 
nonscattering  solution  is  a  function  of  lifetime  (and 
also  fluorophore  concentration),  one  would  expect 
that  in  the  presence  of  scatter,  similar  trends  would 
be  observed.  Figure  5  illustrates  the  values  of  the 
fluorescent  phase  shift  at  each  detector  as  the  life¬ 
time  of  the  heterogeneity  is  changed  from  0.1  ns  to 
1000  ns.  At  the  given  modulation  frequency  (150 


MHz  in  this  calculation),  0,„  first  increases,  reaches  a 
maximum,  and  then  subsequently  decreases  as  t  is 
increased  from  0.1  ns  to  1000  ns.  There  is  a  complex 
interplay  of  the  signal  coming  from  the  object  for 
which  we  are  changing  the  lifetime  and  the  signal 
from  the  background  that  leads  to  the  behavior  de¬ 
scribed.  In  summary,  both  M„,  and  0„,  at  each  de¬ 
tector  are  directly  and  strongly  influenced  by  the 
value  of  lifetime  in  the  heterogeneity. 

Detailed  analytical  expressions  for  fluorescence 
amplitude  and  phase  shift  for  infinite  and  semi- 
infinite  media  with  spherical  heterogeneities  have 
been  provided  by  Li  et  aZ.29  Our  numerical  compu¬ 
tational  methods  provide  a  simulation  of  finite  media 
and  arbitrary  shaped  hidden  objects.  Our  numeri¬ 
cal  results  agree  with  the  general  predictions  pro¬ 
vided  by  Li  et  al ,29 

The  solution  of  the  forward  problem  was  used  as 
inputs  to  the  solution  of  the  inverse  problem  de- 


Fig.  5.  Plot  of  fluorescent  light  phase  shift  at  150  MHz  at  various 
detectors  for  a  heterogeneous  tissue  phantom  with  a  30-mm- 
diameter  object  located  at  the  center  of  the  phantom  for. various 
object  t  values.  Background  lifetime  t  is  1  ns  and  absorption 
coefficients  Tip.0i  n  are  1  x  10  s  mm  '  for  the  background  and  1  x 
10  mm  1  for  the  object. 
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Table  2.  Optical  Properties  Used  to  Generate  the  Simulated  Experimental  Data  as  Inputs  to  the  Inverse  Image-Reconstruction  Algorithm 


Case 

(mm ‘  l) 

M'r, 

(mm1) 

5.A 

0.0 

0.0 

5.B 

1.0  x  10 

1.0  X  10  ; 

5.C 

0.0 

0.0 

5.D 

1.0  X  1(T:' 

1.0  x  10  : 

T 

crM  Gaussian  Noise 

o-(i  Gaussian  Noise 

M-Ni  '  or  |xs  ' 

Background 

Background 

in  Log  of  ac 

in  Phase 

(mm  l) 

(ns) 

(mm'1) 

Amplitude 

(deg) 

1.0 

10.0 

1.0  x  10~5 

0.01 

0.1 

1.0 

10.0 

1.0  x  10~5 

0.01 

0.1 

1.0 

10.0 

1.0  x  10“5 

0.01 

1.0 

1.0 

10.0 

1.0  x  10'3 

0.01 

0.1 

scribed  in  Section  4.  Forward  solutions  in  two  di¬ 
mensions  were  computed  for  the  three  cases  outlined 
in  Tables  2  and  3  for  the  source  and  detector  geom¬ 
etry  in  a  100-mm-diameter  simulated  tissue  phantom 
(Fig.  1).  Values  of  0m  and  Mm  were  computed  at 
each  detector  in  response  to  four  sources  of  excitation 
light  located  at  the  periphery  modulated  at  a  single 
frequency  of  150  MHz.  Consequently,  the  forward 
solution  predicted  80  values  of  0m  and  Mm,  each  at 
detector  i  (i  =  1,  20)  in  response  to  source  k  (k  =  1,  4). 
Gaussian  noise  with  a  standard  deviation  of  0.1°'  (or 
a  liberal  1°)  in  9m  and  1%  in  Mm  was  superimposed  on 
the  solution  of  the  forward  problem.  These  obtained 
data  sets  are  used  as  the  simulated  experiments  as 
input  to  the  reconstruction  algorithm  described  be¬ 
low. 

4.  Inverse  Problem  and  the  Solution  Methodology 

We  begin  our  inverse  calculations  from  a  uniform 
starting  guess  of  a  fluorescence  yield,  (-q^  )  ,  and 

lifetime,  (t)7,  at  all  the  grid  points  within  tKe  tissue 
phantom.  With  the  initial  guess  of  an  optical  prop¬ 
erty  map,  one  can  evaluate  the  complex  fluence  at 
the  various  detector  locations.  There  are  two  ap¬ 
proaches  that  one  can  follow:  one  involves  mini¬ 
mizing  the  difference  between  the  complex  number 
that  is  calculated  from  the  ac  magnitude  and  phase, 
and  the  other  uses  the  ac  magnitude  and  phase 
information  itself.  The  two  approaches  are  equiv¬ 
alent,  and  we  have  chosen  the  second  approach  of 
using  the  ac  magnitude  and  phase.  With  the  start¬ 
ing  guess  of  optical  properties,  we  compute  a  pre¬ 
diction  of  phase  shift  (0,„).  and  log  of  ac  amplitude 
(Af,,,),-  at  each  detector  i.  Values  of  (r|p.„  )  and 

(t);  are  iteratively  adjusted  at  each  grid ’point  to 
minimize  the  error  between  the  predicted  and  sim¬ 


ulated  experimental  values  of  (0m)i  and  (AfJ,-  re¬ 
sulting  from  each  source,  k,  k  =  1,  4.  The  choice  of 
appropriate  functions  to  be  minimized  for  our  re¬ 
constructions  is  now  discussed.  Section  3  de¬ 
scribes  at  length  how  changes  in  the  fluorescent 
yield  and  lifetime  affect  the  measured  values  of  ac 
magnitude  and  phase.  Changes  in  fluorescent 
yield,  T)p.a  ^  ,  have  a  direct  and  strong  influence  on 
the  ac  magnitude  and  only  an  indirect  and  weak 
influence  on  the  phase.  Hence  we  have  chosen  to 
adjust  iteratively  the  values  of  ('np.a,_„,)/  in  order  to 
minimize  the  merit  function,  x,*2,  that  depends  on 
the  ac  magnitude  measurements: 


where  aM  is  the  typical  standard  deviation  of  noise  in 
Mm,  taken  to  be  0.01. 

The  value  of  Mm^.  is  the  simulated  experimental 
value  computed  in  Section  3  with  added  Gaussian 
noise.  The  goal  of  the  algorithm  is  to  minimize  x  2 
by  appropriate  updates  of  (ti|xa  After  an  itera¬ 
tion  of  the  (TiP-ai_m),  update,  values  of  (t)7-  were  ad¬ 
justed  in  the  next  iteration  in  order  to  minimize  a 
second  merit  function,  xT2: 


,-*M2  + 

<T.W  /  \  0 r„  /  ’ 

(6) 

where  ct„  is  the  typical  standard  deviation  of  noise 
in  0,  taken  to  be  1°.  The  values  of  0  are  the 
simulated  experimental  values  computed  m  Section 
3  with  added  Gaussian  noise.  The  above  merit 
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Table  3.  Location  and  Area  of  the  Simulated  Heterogeneities:  Comparison  of  Expected  and  Reconstructed  Values 


Object  1 

Object  2 

Case 

Area 

(mm2) 

Location  (x,  y) 

(mm) 

Area 

(mm2) 

Location  (jc,  y) 

(mm) 

5. A 

5.B 

5.C 

5.D 

706.0  (expected) 
742  (obtained) 
706.0  (expected) 
683  (obtained) 
314.1  (expected) 
381  (obtained) 
706.0  (expected) 
693  (obtained) 

(60,  60)  (expected) 

(61,  59)  (obtained) 

(60,  60)  (expected) 

(59,  58)  (obtained) 

(32.3,  67.7)  (expected) 

(34,  68)  (obtained) 

(60,  60)  (expected) 

(61,  57)  (obtained) 

not  applicable 

not  applicable 

314.1  (expected) 
342  (obtained) 
not  applicable 

not  applicable 

not  applicable 

(67.7,  32.3)  (expected) 

(65,  35)  (obtained) 
not  applicable 
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function  is  chosen  because  changes  in  lifetime  7 
have  a  direct  and  strong  influence  on  both  the  ac 
magnitude  and  phase.  Hence  we  have  chosen  to 
adjust  iteratively  the  values  of  (T).  in  order  to  min¬ 
imize  the  merit  function  that  consists  of  both  ac 
magnitude  and  phase  measurements.  The  goal  of 
this  part  of  the  algorithm  is  to  minimize  x  2  by 
appropriate  choice  of  (t )j  as  described  below!  In 
Eq.  (5),  the  choice  of  the  value  of  < tm  does  not  affect 
the  minimization- update  process.  We  have  cho¬ 
sen  the  typical  noise  values  for  M  and  0,  aM  and  ct 
respectively,  to  be  factors  for  converting  to  dimem 
sionless  numbers.  In  fact  crw  and  o-„  are  simply 
scale  factors  that  assess  confidence  in  measure¬ 
ment.  For  example,  when  crH  is  smaller  for  a  given 
(tm,  the  0  residuals  are  weighted  more  heavily  in  the 
inversion  as  compared  with  the  M  residuals. 
Thus,  factor  a  is  necessary  to  normalize  the  indi¬ 
vidual  sum  of  residuals  of  M  and  0  so  the  two  can  be 
added  in  Eq.  (6). 

To  update  values  of  (rip.a . )•  and  (t)  ,  one  needs  the 

Jacobian  matrices  that  describe  the  sensitivity  of  the 
detector  response  at  position  i  to  changes  in  (t)  and 
at  each  grid  point,  j.  The  elements  of  the 
three  Jacobian,  matrices  employed,  j(M  rm  ) 
m,  t),  and  =7(0,  t),  are  given  byL’-  W/ 
. —  (dMj/dTj),  and j\j  =  (<50;/Jt  ),  respec¬ 
tively.  One  calculates  these  elements  by  solving  the 
forward  problem  four  times  for  each  grid  point,  j  to 
obtain  Mm  i  and  0m  •  calculated  with  (t)-  and  (t  4-  At) 
and  with  and  +  Atip^J;.  From 

the  least-squares  minimization,  one  can  show  that 
the  update  in  Tip.a  ^  and  t  can  be  calculated  by  using 
an  algorithm  (Newton’s  method)  similar  to  that  sug¬ 
gested  by  Yorkey  et  al.30  for  the  reconstruction  of 
images  obtained  by  electrical  impedance  tomogra¬ 
phy. 

Equations  (7)  and  (8)  provide  updates 
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and  [At]  to  the  yield  [=fpa  _  ]  and  lifetime 
LxJ  vectors  at  each  iteration.  The  upcfate  of  Op  ] 
is  based  on  the  ac  amplitude  data,  whereas  that  of7?l 
is  based  on  both  the  ac  amplitude  and  phase  data. 
This  follows  from  the  discussion  in  the  previous  sec¬ 
tion  on  the  forward  problem.  Mm^  and  Mm  are  the 

experimentally  simulated  and  calculated  vectors 
consisting  of  the  log  of  the  ac  amplitude  at  each 
of  the  i  detectors,  respectively.  Because  of  the  ill- 


conditioned  nature  of  the  Jacobian  matrices,  terms 
"J 1  or  M  are  added  (/  is  an  identity  matrix)  to  make 
the  matrices  more  diagonally  dominant  and  the  so¬ 
lution  of  the  algebraic  equations  more  robust.  Pa¬ 
rameters  or  \2  are  adjusted  by  means  of  a 
Marquardt-Levenberg  type  of  algorithm.31  Lower 
triangular- upper  triangular  decomposition  and 
backsubstitution  are  employed  to  solve  the  simulta¬ 
neous  linear  algebraic  equations  given  in  Eqs.  (7)  and 
(8L31  At  each  iteration,  the  merit  functions,  the  Ja- 
cobian  matrices,  and  the  updates  in  the  fluorescent 
yield  and  lifetime  are  evaluated  and  the  iterations 
are  continued  until  the  convergence  criterion  is  met. 
Convergence  is  achieved  when  any  of  the  following 
three  quantities,  i.e.,  x2,  change  in  v2  in  successive 
iterations,  and  relative  change  in  x2  in  successive 
iterations, ^  is  lower  than  a  predetermined  value  of 
1.0  X  10 

5.  Results  and  Discussion 

The  performance  of  FLI  by  using  the  inversion  algo¬ 
rithms  described  above  is  depicted  in  Figs.  6-9  for  the 
case  studies  listed  in  Tables  2  and  3.  Simulated 
experiment  5. A  was  designed  to  reconstruct  (x)  and 
with  no  absorption  that  is  due  to  nonfluo- 
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(b) 


Fig.  7.  Reconstructed  spatial  map  of  fluorescence  (a)  yield, 
(b)  lifetime,  t,  on  a  two-dimensional,  17  X  17  grid  for  the 
case  described  in  Subsection  5.A.  ^  at  the  excitation  wave- 
length  not  accounting  for  ftuorophore  absorption  is  zero:  at  the 

emission  wavelength  is  also  zero.  The  average  values  of  "-nix 
and  t  within  the  object  were  1  x  10  '  mm  '  and  1  ns,  respectively 
(expected)  and  0.93  x  10  :l  mm  1  and  1.03,  respectively  (recon¬ 
structed).  Spurious  unphysically  high  values  of  and  r 

have  been  replaced  by  the  average  background  fluorescence  yield 
and  lifetime,  respectively,  obtained  from  the  inversion. 


100r  xlCT3 


(b) 

Fig.  8.  Reconstructed  spatial  map  of  fluorescence  (  a)  yield,  ^ 

(b)  lifetime,  t,  on  a  two-dimensional,  17  X  17  grid  for  the  case 
described  in  Subsection  5.B.  ^  at  the  excitation  wavelength  not 
accounting  for  fluorophore  absorption  is  1  X  10“;J  mm”1;  jia  at  the 
emission  wavelength  is  zero.  The  average  values  of  and  r 

within  the  object  were  1  X  10  {  mm  1  and  1  ns,  respectively  (ex¬ 
pected)  and  0.8  x  10  mm  1  and  0.7  ns,  respectively  (reconstruct¬ 
ed).  Spurious  unphysically  high  values  ofij^  and  r  have  been 
replaced  by  the  average  background  fluorescence  yield  and  lifetime, 
respectively,  obtained  from  the  inversion. 


rescing  chromophores,  whereas  simulated  experi¬ 
ment  5.B  included  a  finite  chromophore  absorption  of 
excitation  and  fluorescent  light  that  could  be  consid¬ 
ered  physiological.  Finally,  simulated  experiment 
5.C  evaluated  the  ability  to  determine  the  location  of 
two  hidden  objects  within  the  tissue  phantom 
whereas  5.D  examined  the  reconstruction  when  the 
uptake  ratio  uatio  of  fluorescent  yield  in  heterogene¬ 
ity  to  that  in  the  background)  was  20. 

A.  Single  Heterogeneity  with  Optical  Contrast  from  (x) 

and  (TlAa,.  J,  with  Absorption  of  Excitation  and  Fluorescent 
Light  by  Chromophores 

To  calculate  the  experimental  data  for  the  first  case 
the  fluorescence  yields,  xux„  ,  for  the  background 


and  a  single  object  were  chosen  as  1  X  10~5  mm-1 
and  1  x  10  mm  \  respectively,  and  the  lifetimes, 
t,  tor  the  background  and  the  object  were  chosen  as 
10  and  1  ns,  respectively.  During  the  inverse  im- 
age  reconstruction,  no  a  priori  knowledge  of  either 
the  object  location  or  the  background  fluorescence 
properties  was  assumed  and  a  uniform  guess  of  1  X 
10  mm  1  and  10  ns  was  given  for  -rjp,a  and  x, 
respectively.  Convergence  was  achieved'  in  fewer 
than  50  iterations  (computational  time  on  a  Sun- 
Sparc  10  was  2  h)  for  a  two-dimensional  17  x  17 
grid.  The  average  values  of  iqp.,,  and  x  in  the  grid 
points  that  occupy  the  simulated  object  converge 
within  50  iterations  to  the  values  of  T]|ia  =  0.93  X 
10  mm  and  x  =  1.03  ns,  which  are  close  to  the 
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Fig.  9.  Reconstructed  spatial  map  of  fluorescence  yield,  pp 
on  a  two-dimensional,  33  x  33  grid  for  the  case  described  in  Sub¬ 
section  5.C.  p.a  _  at  the  excitation  wavelength  not  accounting  for 
fluorophore  absorption  is  zero;  p.^  at  the  emission  wavelength  is 
also  zero.  The  Gaussian  noise  that  was  introduced  in  the  phase 
had  a  standard  deviation  of  1  °.  The  object  locations  and  sizes  are 
recovered  correctly.  The  average  values  of  pp  within  the  two 
objects  were  top  left,  1  x  10~3  mm  1  (expected)  and  2  x  10  3 
mm'1  (reconstructed);  bottom  right,  2  x  10  '3  mm1  (expected) 
and  1.8  x  10  mm  (reconstructed).  Spurious  unphysically 
high  values  of  pp„>_.,ji  have  been  replaced  by  the  average  back¬ 
ground  fluorescence  yield  obtained  from  the  inversion. 


correct  values  of  tipt0(  „  =  1  X  1(T3  mm-'  and  t  = 
1  ns,  as  shown  in  Figs.  6(a)  and  6(b)  and  listed  in 
Table  4.  Figures  7(a)  and  7(b)  illustrate  the  recon¬ 
structed  images  of  T|p,0^  and  t,  respectively,  and 
are  representative  of  the  expected  images.  The 
images  were  smoothed  by  interpolation  in  this  and 
subsequent  cases.  During  the  reconstruction,  it 
was  observed  that  some  grid  points,  often  on’  or 
close  to  the  periphery,  had  unphysically  high  values 
of  the  yield  as  well  as  lifetime.  Typically,  these 
grid  points  were  loners  and  were  surrounded  by 
grid  points  with  reasonable  yield  and  lifetime  val¬ 
ues.  The  high  values  of  lifetime  lead  to  a  lower 


Table  4.  Fluorescence  Lifetime  t  and  Yield  ppa,  of  the  Simulated 
Heterogeneities:  Comparison  of  Expected  and  Reconstructed  Values 


Case 

"nM'r/j  (Object) 

(mm  J) 

t  (Object) 

(ns) 

5.A 

1.0  X  10“ 3  (expected) 

1.0  (expected) 

5.B 

0.93  x  10  “■*  (obtained) 

1.03  (obtained) 

1.0  x  10  (expected) 

1.0  (expected) 

5.C 

2.1  X  10  3  (obtained) 

4.1  (obtained) 

Top  left  object 

1.0  X  10' 3  (expected) 

1.0  (expected) 

Bottom  right 
object 

2  x  10- 3  (obtained) 

4.1  (obtained) 

2.0  x  10  3  (expected) 

2.0  (expected) 

5.D 

1.8  x  10  3  (obtained) 

3.5  (obtained) 

2.0  X  10  4  (expected) 

1.0  (expected) 

7.0  X  10  4  (obtained) 

4  (obtained) 
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magnitude  of  fluorescence  generation  and  offset  the 
effect  of  high  yield  values.  We  believe  that  addi¬ 
tional  criteria,  such  as  smoothness  of  the  recon- 
stiucted  maps,  may  alleviate  the  spurious 
background  peak  values.  In  this  studv,  these  spu¬ 
rious  values  were  replaced  by  the  average  back¬ 
ground  fluorescence  yield  and  lifetime  obtained 
from  the  inversion.  In  the  subsequent  reconstruc¬ 
tions,  spurious  values  were  similarly  handled. 

The  average  values  of  T)p.n .  in  the  grid  points 

that  occupy  the  simulated  background  converge 
within  50  iterations  to  9  X  10“5  mm'1  in  compar¬ 
ison  with  the  correct  value  of  1  x  10_r>  mm' 1  (data 
not  shown  for  brevity).  The  value  of  background  r 
converges  to  5.4  ns,  which  is  also  far  from  the  cor¬ 
rect  value  of  10  ns.  Both  discrepancies  are  attrib¬ 
uted  to  the  fact  that  most  of  the  signal  contribution 
is  due  to  the  object,  with  the  object-to-background 
uptake  ratio  being  100:1.  The  dependence  of  the 
final  images  on  the  choice  of  the  initial  guess  was 
examined  by  providing  an  initial  uniform  guess  of 
1  X  10  mm  and  10  ns  for  T]p.a  and  t,  respec¬ 
tively.  This  resulted  in  images  similar  to  those 
obtained  with  the  first  guess. 

The  location  of  the  heterogeneity  was  identified  as 
consisting  of  all  the  grid  points  with  T)|xa  higher 
than  35%  (arbitrarily  chosen)  of  the  peak  value  of 

. [Fig-  7(a)],  The  average  of  the  coordinates  of 

all  the  identified  object  grid  points  was  the  position 
(60.8,  58.5),  which  is  close  to  position  (60,  60)  that  was 
used  to  simulate  the  experimental  data.  As  listed  in 
Table  3,  the  area  of  the  heterogeneity  based  on  our 
arbitrary  definition  for  identification  was  742  mm2, 
close  to  that  used  to  generate  our  simulated  experi¬ 
mental  data.  In  addition,  upon  inspection  of  Fig. 
7(a),  one  can  note  that  values  of  _  for  the  object 
on  grid  points  closer  to  the  center'" of  the  phan¬ 
tom  were  higher  than  those  toward  the  periphery. 
This  is  due  to  a  reduction  of  signal  contributions  of  grid 
points  located  farthest  away  from  source  and  detec¬ 
tors.  This  is  reflected  in  a  smaller  value  in  the 
Jacobian  matrices  and  results  in  a  poorer  reconstruc¬ 
tion  at  the  center  of  the  phantom. 

B.  Single  Heterogeneity  with  Optical  Contrast  from  (t). 
and  with  Absorption  of  Excitation  and  Fluorescent 

Light  by  Chromophores 

The  same  hidden  object  location  as  well  as  optical 
parameters  were  used  as  described  in  Subsection 
5. A,  except  that  a  uniform  background  chro- 
mophore  absorption  coefficient  of  the  excitation 
light,  p.Qr  of  1  X  10  3  mm  1 ,  and  of  the  fluorescent 
light,  p,  of  1  X  10“3  mm-1,  was  used  to  generate 
the  simulated  experimental  data.  Although  exci¬ 
tation  light  propagation  was  not  employed  for  im¬ 
age  reconstruction,  we  considered  this  optical 
property  known  to  estimate  the  best  possible  per¬ 
formance  for  inverse  image  reconstruction  under 
physiological  conditions.  The  two-dimensional  re¬ 
constructed  spatial  map  of  the  fluorescence  yield, 
Wo,.™'  an(3  lifetime,  t,  are  shown  m  Figs.  8(a)  and’ 

8(b),  respectively.  The  image  quality  is  slightly 


lltZ°rfW^h  ue^ect  *?  both  the  size  as  well  as  the 
shape  of.  the  hidden  object.  As  shown  in  Table  3 

the  mean  value  of  location  of  the  object  according  to 

criterion  based  on-rm  .  occurred  at  position 

(59  58),  consistent  with  the  conditions  used  to  sim¬ 
ulate  the  experimental  data.  The  dimension  of  the 

°?  °Ur  arbitrary  definition  for 
identification  (all  grid  points  with  ^  higher 

than  35%  of  the  maximum)  was  683  mm5',’ which  is 
close  to  that  used  to  generate  our  simulated  exper¬ 
imental  data.  The  average  values  of  t|(x0  and  t 

m  the  grid  points  that  occupy  the  simulated  object 
converge  within  50  iterations  to  the  values  of  t„x  = 

bip-hpr  th  fTm  i  3nC*  T  =  4  ns>  which  are  slightly 
higher  than  the  values  used  to  generate  the  simulated 

experimental  data  (see  Table  4).  The  average  values 

l  Tj‘  r  aPd  T  ln  the  grid  points  that  occupy  the  sim¬ 
ulated  background  converge  within  50  iterations  to 
values  similar  to  those  reported  in  Case  5.A. 

C.  Two  Heterogeneities  with  Optical  Contrast  from  (T). 
and  (ruia,_J/  with  No  Absorption  of  Excitation  and  ' 
Fluorescent  Light  by  Chromophores 

In  Case  5.C,  the  same  optical  parameters  were  used 
in  forward  calculations  as  described  in  Subsection 
5.A,  except  that  the  fluorescence  yields,  pixa  ,  for 
objects  1  and  2_were  chosen  as  1  x  1<T3  mm'1  and 
1*  }°  .  mm  -  respectively,  and  lifetimes,  t,  for 
the  objects  were  chosen  as  1  ns  and  2  ns  resnec- 
t,,Ve  y;,  ,Two  20-mm-diameter  circular  objects  were 
Seda,0nfKa  dla®°ilaI  at  the  coordinates  shown  in 
phantom Wlthm  3  100'mm'diameter  circular  tissue 

Again  during  the  reconstruction,  no  a  priori  knowl¬ 
edge  of  either  the  object  location  or  the  background 
fluorescence  properties  was  assumed  and  the  values 
of  yield  nP-a,-„, and  lifetime  t  were  found  at  all  the  grid 
points  on  a  33  X  33  grid.  The  two-dimensional  re¬ 
constructed  spatial  map  of  the  fluorescence  yield 
is  shown  in  Fig.  9.  The  object  locations  (*,  y) 
obtained  are  given  in  Table  3  and  match  well  with  the 

tTdatanS  Thf t0  generaS the  simulated  experimen- 
tal  data.  The  areas  of  the  objects  from  the  recon- 

structed  image  (ah  grid  points  with  vtia  higher 
than  35%  of  the  maximum)  were  381  mm2  (top 
left  object  1)  and  342  mm2  (bottom  right,  object  2) 
slightly  larger  than  inputs  to  the  forward  problm 
The  average  values  of  ^  and  t  in  the  grid  poinTs 

a  occupy  the  simulated  object  converge  within  100 
iterations  to  the  values  of  T|pa  =  2  X  10“3  mm"1 


40  60 

x  coordinate  (mm) 

Fig.  10.  Reconstructed  spatial  map  of  fluorescence  yield  Tm 

senctaionW5  DlmenSi0ni’;.?3  X  33  ^  f°r  the  C3Se  desCribed  in  Sub-’ 
fluorophore  n.aeco  ting for 

wavelength  is  the  same.  The  average  value 


mm  1  (re- 


-  -  - —  avcidge  value  OI  711X„ 

o  ject  was  2  X  10'  1  mm"1  (expected)  and  7.0  X  l(p  mr 
constructed).  Spurious  unphysically  high  values  of  have 

been  replaced  by  the  average  background  fluorescence  jrteld  ob¬ 
tained  from  the  inversion.  ^ 


10  5  mm  \  whereas  the  value  of  the  background  life¬ 
time  was  poorly  reconstructed. 

D  Single  Heterogeneity  with  Optical  Contrast  from  (t) 

iTnht  rW'  Wlth  Absorption  of  Excitation  and  Fluorescent 
Light  by  Chromophores  and  an  Uptake  Ratio  of  20 

The  same  hidden  object  location  as  well  as  optical 

excTnTthat  I®  d?scribed  in  Subsection  5.A, 

except  that  a  uniform  background  chromophore  ab¬ 
sorption  coefficient,  of  1  X  1(T3  mm"1  as  well  as 
P-a„,  at  the  emission  wavelength  of  1  x  10-3  mm-1 

daeta  anddr°  g6nerate  ■the  simulated  experimental 
data,  and  Gaussian  noise  was  added  as  discussed 

a£dVtb  VS°; the  values  of  for  the  background 
aad4the  oble,cfc  'yere  chosen  as  1  X  10~5  mm"1  and  2  X 

1 V  mm  trine  rninncr  _ 1  1  ,  .  _ 


and  ^ .  ~  L8  x  10  3  mm  1  for  objects  1 

and  2,  respectively  which  again  is  close  to  the  expected 
values  of  1  X  1(T3  mm"1  and  2  X  1(T3  mm-i  Life, 
times  of  the  two  objects  were  found  to  be  4.1  and  3.5  ns 
respectively,  whereas  the  expected  values  were  1  and  2 
ns  respectively.  The  quantitative  values  of  m 
and  t  obtained  from  the  inverse  solution  are  currently 
unsatisfactory,  and  research  is  in  progress  to  improve 
the  solution  procedure.  The  background  vZe  of 
agrees  well  with  the  expected  value  of  1.0  X 


If)-4  m™-l  4.U  •  •  ■  '  *“  **“■“  auu  &  A 

Apoin^uu  ’  tuUS  glvlng  an  uptake  ratio  of  20. 
Pmni  ’  lth°ugh  excitation  light  propagation  was  not 
employed  for  image  reconstruction,  we  considered 

thf  KPttCaI  Pr??erty  t0  be  known  in  order  to  estimate 
the  best  possible  performance  for  inverse  image  re- 

constructmn  under  physiological  conditions. &  The 
two-dimensional  reconstructed  spatial  map  of  the  flu¬ 
orescence  yield,  Tiitaw  is  shown  in  Fig  10  The 
rnicige  quality  with  respect  to  both  the  size  as  well  as 

Tableh3P tb°f  ^  b^Z  °bject  is  g00cL  As  shown  in 
Table  3,  the  mean  value  of  location  of  the  object  ac- 

posS  (eiTyf1 teri°n,baf ' 1  °un  occurred  at 

to  simnlit  Vk?)’  consistent  with  the  conditions  used 
of  e  the  exper,lmentaI  data.  The  dimension 

of  the  heterogeneity  based  on  our  arbitrary  defini- 
tion  for  identification  (all  grid  points  with  ¥ 
lgher  than  35%  of  the  maximum)  was  693  mm2" 
which  is  close  to  that  used  to  generate  our  simu¬ 
lated  experimental  data.  The  average  values  of 
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an<*  T  in  the  grid  points  that  occupy  the  sim¬ 
ulated  object  converge  within  50  iterations  to  the 
values  of  =  7  x  1(T4  mm'1  and  t  =  4  ns 

which  are  slightly  higher  than  the  values  used  to 
generate  the  simulated  experimental  data  (see  Table 
4).  The  average  values  of  iip*  and  t  in  the  grid 
points  that  occupy  the  simulated  background  con¬ 
verge  within  50  iterations  to  values  similar  to  those 
reported  in  Case  5.A.  Nonetheless,  it  is  shown  that 
image  reconstruction  may  be  successful  even  in  the 
case  of  uptake  ratio  as  small  as  20. 

The  performance  of  the  inverse  imaging  algorithm 
may  be  improved  by  including  additional  information 
obtained  from  multifrequency  measurements  and 
from  excitation  wavelength  measurements  In  ad¬ 
dition,  the  appropriate  weighting  of  grid  point  contri¬ 
butions  based  on  the  signal  strength  may  also 
improve  reconstruction  results. 

6.  Conclusions 

The  fluorescent  yield  and  lifetime  of  an  endogenous 
fluorophore  may  be  sensitive  to  local  environments 
providing  specificity  for  contrast  of  diseased  over 
normal  tissues  and  optimal  detection  of  disease  by 
using  optical  techniques.32  Our  simulated  experi¬ 
ments  show  that  it  is  possible  to  reconstruct  the 
fluorescent  yield  and  lifetime  of  embedded  fluoro- 
phores  in  tissue-mimicking  scattering  media  from 
frequency-domain  measurements  of  fluorescent 
phase  shift  and  ac  amplitude  or  amplitude  modula¬ 
tion.  As  we  have  shown  in  Subsection  5.C,  the 
reconstruction  of  lifetime  can  be  problematic  by  using 
the  reemitted  fluorescent  signal  at  one  modulation 
frequency.  This  is  in  agreement  with  O’Leary  et 
at.,33  who  resorted  to  the  need  for  measurements  in 
the  absence  of  fluorophore  (background  properties)  as 
well  as  in  its  presence  for  lifetime  reconstruction  as 
suggested  from  multipixel  measurements  for  photon 
migration  imaging.34  Currently,  we  are  experimen¬ 
tally  investigating  the  implementation  of  FLI  by  us¬ 
ing  both  excitation  and  fluorescent  wavelengths  as 
well  as  multifrequency  measurements  to  improve  in¬ 
verse  solutions  that  predict  both  fluorescent  yield  and 
lifetime. 

Appendix  A:  Nomenclature 

c,  velocity  of  light 
D(r)}  optical  diffusion  coefficient 
f >  modulation  frequency 
/,  identity  matrix 

i,  detector  number,  i  =  1,  20 
J,  Jacobian  matrix  relating  the  sensitivity  of  op¬ 
tical  parameters  to  the  detector  response 

j,  grid  point  number 

jij,  individual  elements  of  Jacobian  matrix  J 

k,  source  number,  k  =  1,4 

M,  log  of  ac  amplitude  of  modulated  fluorescent 
light 

r ,  position 

Sir,  a,),  source  term  for  the  modulated  light  at  position 
r  and  frequency  u> 

2270  APPLIED  OPTICS  /  Vol.  36,  No.  10/1  April  1997 


Greek 

2 

X  ,  merit  function  representing  the  least- 
squares  error 

^(r,  w),  complex  number  representing  the  photon 
flux  in  the  frequency  domain  at  position  r 
and  frequency  co 

ti,  quantum  yield  of  the  fluorescent  probe  or 
dye 

p.*,  average  absorption  coefficient  of  the  tissue 
PQ/n,  absorption  coefficient  of  the  fluorescence 
light  by  both  the  nonfluorescing  chro- 
mophores  and  fluorophore 
Par,  absorption  coefficient  of  the  excitation  light 
by  both  the  nonfluorescing  chromophores 
and  fluorophore 

^arn_i  absorption  coefficient  of  the  fluorescence 
light  by  both  the  nonfluorescinff  chro¬ 
mophores 

absorption  coefficient  of  the  excitation  light 
by  the  nonfluorescing  chromophores 
absorption  coefficient  of  the  excitation  light 
by  the  fluorophore 

Ps  >  effective  scattering  coefficient  of  the  tissue 
0,  phase  shift  of  the  modulated  light  wave 
with  respect  to  the  modulated  wave  at  the 
source 

o',  standard  deviation  of  the  Gaussian  noise 
representing  the  experimental  uncertainty 
T(r),  lifetime  of  the  activated  probe  or  dye  at 
location  r 

co,  angular  modulation  frequency,  given  by 
2irf 

Subscripts 

obs,  observed  or  experimental  data 
x,  excitation  light 
m,  fluorescence  or  emission  light 

Appendix  B:  Jacobian  Matrices 

We  provide  more  details  about  the  elements  of  the 
Jacobian  matrices  introduced  in  Section  4.  Element 
Jij  of  a  Jacobian  matrix  represents  the  sensitivity  of 
the  response  of  detector  i  to  changes  in  the  optical 
property  at  grid  pointy.  Here  as  an  example  we  show 
the  response  of  detector  16  for  source  A  (see  Fig  1)  to 
changes  in  t  and  ^  at  all  the  grid  points  for  the 
case  described  in  Subsection  5.A  for  the  first  iteration 
during  inversion  when  |xa  _  was  chosen  as  zero.  Sim¬ 
ilar  results  are  observed' for  other  cases.  Figures 
11(a)  and  ll(b)_show  the  elements  of  Jacobian  matri- 
ces  J{9  x)  and  J(M,  t),  respectively,  where  t  at  each  of 
the  grid  points  was  increased  by  5%.  Similarly,  Fig. 

11(c)  shows  the  elements  of  Jacobian  matrix  ’j{M, 
where  at  each  of  the  grid  points  was 

increased  by  1%.  Most  of  the  elements  of  j(6,  t) 
and  are  positive,  whereas  most  of  the  ele¬ 

ments  of  J(M,  t)  are  negative.  This  is  according  to  our 
expectation  of  systems  with  no  scattering.  The  de¬ 
pendence  on  scatter  has,  of  course,  been  taken  into 
account  in  the  above  example. 


0.15 


y-coordinate  (mm)  0  0 

x-coordinate  (mm) 


(a) 


x  10"4 


y-coordinate  (mm)  0  0 

x-coordinate  (mm) 


(C) 

Fig.  11.  Jacobian  fa)  J(fl,  t),  (b)  j(M,  t),  (c)  J(M,  r,p.a  )  for  the 
case  described  in  Subsection  5.A  for  source  A,  detector1 16,  and 
iteration  1.  During  the  computation  of  the  Jacobians,  the  vklues 
oft  for  (a)  and  (b)  and  ^  for  (c)  at  each  individual  grid  point 
were  increased  by  5%,  5 9c,  and  1%,  respectively. 
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Role  of  higher-order  scattering 
in  solutions  to  the  forward  and  inverse 
optical-imaging  problems  in  random  media 
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biomedical  optical  imaging. 


1.  Introduction 

With  the  development  of  near-infrared  (NIR)- 
emitting  laser  diodes  and  the  realization  that  NIK 
light  can  travel  several  centimeters  through  tissue, 
numerous  groups  have  embarked  upon  development 
of  optical  imaging  as  a  new  medical  imaging  modal¬ 
ity.  Although  approaches  such  as  monitoring  the 
vanishingly  small  component  of  coherent  light  with 
sophisticated  techniques  of  Kerr  filtering,1  time  gat¬ 
ing,2  and  conservation  of  light  polarization3  vary, 
other  approaches  focus  on  monitoring  the  predomi¬ 
nant  optical  signal  re-emitted  from  tissues,  the 
multiply  scattered  signal.  Continuous  wave,4  time 
domain,5  and  frequency  domain6-8  measurements  of 
multiply  scattered  light  have  been  performed  in  sim¬ 
ulated  and  experimental  tissue  phantom  studies  as 
well  as  in  human  studies.9  However,  the  use  of 
these  measurements  for  reconstructing  maps  of  in¬ 
ternal  tissue  optical  properties  for  diagnostic  imaging 
has  been  problematic  because  the  geometric  correla¬ 
tion  between  incident  and  detected  radiation  is  de¬ 
stroyed  by  multiple  scattering. 

While  some  investigators  have  used  direct  image 
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reconstruction  approaches  that  use  measured  optical 
data  to  directly  form  an  image,9-11  others  have 
sought  full  solution  to  the  inverse-imaging  problem, 
which  relates  external  time-dependent  measure¬ 
ments  made  at  the  periphery  to  an  optical  property 
map  of  interior  volumes  through  the  optical  diffusion 
equation.6-8-12-16  Two  approaches  for  solving  the  in¬ 
verse  solution  have  been  adopted.  In  the  first  ap¬ 
proach,  the  fluence  associated  with  a  propagating 
photon  density  wave  launched  at  source  position  ps 
and  detected  at  the  tissue  periphery  at  position  pd 
[denoted  <fr(p„  pd)]  is  related  to  the  fluence  assumed 
in  the  absence  of  any  optical  heterogeneities  [denoted 
cj>  (ps,  pd)]  and  the  internal  optical  property  pertur¬ 
bation  map  A(p)  through  a  Fredholm  integral  equa¬ 
tion  of  the  first  kind: 


<t>(p„  pd)  =  <tW(ps-  Pi)  + 


G(p,  Pd)^o(Ps>  p)A(p)3  p,  (I) 


where  G(p,  pd)  is  the  Green’s  function  to  the  diffusion 
equation, 


Gip,  Pd )  =  T~r~ - ^  exp{i[(-pi>  +  iw/cn)/D]1/2(p  Pd) I. 

describing  the  propagation  of  light  from  position  p  to 
the  detector  at  pd,  and  <t>0  is  the  incident  wave  im¬ 
pinging  at  position  p.  If  one  assumes  that  the  inci¬ 
dent  wave  impinging  upon  the  position  p  can  be 
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approximated  by  the  fluence  predicted  in  the  absence 
of  heterogeneities  (i.e.,  $0  — *  ^ino  known  as  the  Bom 
approximation),  upon  measurement  of  fluence  at  a 
variety  of  source- detector  separations  and  upon  dis¬ 
cretization  of  Eq.  (1),  an  optical  map  of  perturbations 
A(p)  can  be  obtained  to  a  first-order  approximation.15 
However,  the  Bom  approximation  O0  — ►  <J>inc  does  not 
account  for  second-  and  higher-order  effects  that 
might  arise  from  the  rescattering  of  photon  density 
waves  associated  with  neighboring  inhomogeneities. 
In  addition,  this  approach  assumes  that  perturba¬ 
tions  in  optical  properties  at  a  position  p >-  do  not 
influence  the  propagation  of  light  from  position  p,+1 
to  detector  pd  as  described  by  the  Green’s  function 
G(p,  ptf). 

Investigators  use  the  Bom  iterative  method  (BIM) 
to  account  for  strong  perturbations  and  for  second-  and 
higher-order  scattering  effects  by  using  the  first-order 
approximation  of  optical  property  perturbations  A(p)  to 
iteratively  calculate  4>0  or  the  fluence  incident  at  po¬ 
sition  p.  Yao  et  al. 16  have  shown  that  the  BIM  tends 
to  compensate  for  the  underprediction  of  the  scat¬ 
tering  and  absorption  properties  of  a  single  heter¬ 
ogeneity  in  an  otherwise  uniform  medium  that 
would  occur  when  the  Bom  approximation  (or  single 
iteration)  is  used.  In  addition  to  the  BIM,  the  dis¬ 
torted  Bom  iterative  method  (DBIM)  represents  a 
refinement  in  that  it  also  recompiles  the  Green’s 
functions  G(p,  pd)  to  reflect  changes  in  light  propaga- 
,  tion  and  it  should  speed  convergence.  Yet  to  date 
there  has  been  no  investigation  that  addresses  how 
these  iterative  reconstruction  algorithms  affect  the 
resolution  or,  more  precisely,  how  these  linear  per¬ 
turbation  approaches  can  be  used  to  accurately  image 
closely  positioned,  multiple  heterogeneities  between 
which  nonlinear  second-  and  higher-order  scattering 
effects  exist. 

In  contrast  with  this  inversion  approach,  Jiang  et 
aL 7-8  and  Arridge  et  aL14  have  used  full  numerical  so¬ 
lution  to  the  diffusion  equation  that  describes  the  in¬ 
terdependence  of  voxel  optical  properties  and  their 
contribution  to  the  re-emitted  optical  signal  detected 
from  cw  and  frequency  domain  measurements.  With 
this  approach,  the  second-  and  higher-order  scattering 
arising  from  multiple  inhomogeneities,  which  are  not 
accounted  for  in  the  first  iterations  of  Bom  and  Rytov 
approximations,  are  incorporated  in  the  forward  and 
inverse  solutions.  The  inversion  consists  of  the  relat¬ 
ing  of  spatial  changes  of  optical  properties  on  detected 
frequency  domain  measurements  through  numerical 
solution  of  the  diffusion  equation  and  then  solving  an 
update  to  the  map  of  optical  properties  from  the  dif¬ 
ference  between  the  signals  that  are  measured  and 
those  that  are  predicted  by  the  forward  solution  to  the 
diffusion  equation.  Note  that  convergence  on  the  op¬ 
tical  properties  is  achieved  with  these  numerical  ap¬ 
proaches  and  not  with  the  iterative  Bom  and  Rytov 
approaches.  Nonetheless,  while  the  computational 
investment  of  iterative  but  analytically  based  inver¬ 
sions  is  less  than  that  of  full  numerical  solution  of  the 
diffusion  equation,  the  relative  performance  of  these 


two  inverse-solution  approaches  has  yet  to  be  evalu¬ 
ated. 

For  this  reason  we  embarked  on  a  study  to  assess 
the  contributions  of  neighboring  heterogeneities  by 
using  experimental,  numerical,  and  analytical  com¬ 
putations  of  scattered  photon  density  waves  from 
perfect  light-absorbing  cylinders.  Specifically  we  ex¬ 
perimentally  and  computationally  monitor  the  multi¬ 
ple  scattering  of  photon  density  waves  between  two 
neighboring  perfect  cylindrical  absorbers  embedded  in 
a  tissue-mimicking  scattering  medium  to  assess  the 
higher-order  perturbation  effects  on  a  re-emitted,  de¬ 
tected  photon  density  wave.  In  the  following  we 
briefly  review  the  theory  of  higher-order  perturbation 
effects  as  predicted  from  the  experimental  frequency 
domain  measurements  as  well  as  from  the  Helmholtz 
equations.  In  addition,  we  present  experimental 
measurements  and  numerical  solutions  of  the  optical 
diffusion  equation  that  show  the  errors,  which  can  be 
significant,  introduced  by  neglecting  second-order  ef¬ 
fects.  These  errors  can  define  the  limits  of  resolu¬ 
tion  for  two  perfect  absorbers  in  the  inverse-solution 
approaches  that  do  not  use  BIM,  DBIM,  or  the  full 
numerical  solution  of  the  optical  diffusion  equation. 
We  comment  on  these  effects  when  contrast  is  caused 
by  mechanisms  other  than  a  perfect  absorber. 

2.  Theoretical  Background  for  Perturbations 
Associated  with  Diffuse  Photon  Density  Waves 

Our  work  to  assess  the  contributions  of  higher-order 
perturbations  has  been  motivated  by  the  work  of 
Boas  et  al.15  Their  approach  for  image  reconstruc¬ 
tion  from  diffusely  propagating  photon  density  waves 
uses  analytical  expressions  to  describe  the  complex 
fluence  of  a  propagating  photon  density  wave  <fc(p)  in 
the  presence  of  m  heterogeneities  by  superposition  of 
incident  <l>inc(p)  and  scattered  <J>scatt*n(p)  waves,  from 
body  k  and  of  the  order  of  n: 

*(P)  =  <*>i„c(p)  +  2  2  4W*"(p).  (2) 

n= 1 *=1 

The  second-order  effects  caused  by  the  scattering  of  a 
photon  density  wave  between  two  objects  are  illus¬ 
trated  in  Fig.  1,  from  the  work  of  Boas  and  co-workers, 
and  delineate  the  contribution  of  these  effects  to  the 
detected  photon  density  wave  4>(p).  Second-  and 
higher-order-scattering  effects  arise  from  the  rescat¬ 
tering  of  an  incident  wave  between  the  two  bodies. 
Boas  et  al.  assume  that  second-order-scattering  effects 
(denoted  by  the  dotted  line  in  Fig.  1)  are  negligible 
under  most  circumstances.  We  explore  this  assump¬ 
tion,  using  experimental,  numerical,  and  analytical 
predictions  of  second-order  contributions. 

A.  Analysis  of  Experimental  Measurements  of  Photon 
Density  Waves 

To  measure  second-order  contributions  we  experi¬ 
mentally  measured  the  detected  photon  density  wave 
in  the  presence  of  one  light-absorbing  object  alone, 
object  two  alone,  4>*2;  and  with  both  objects 
present,  as  a  function  of  object  positions  p1  and 
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Fig.  1.  Schematic  illustrating  the  incident  wave  <t>inc  originating 
from  the  source  at  position  p,  (dashed  lines);  the  first-order  scat¬ 
tered  wave  arising  from  the  first  heterogeneity  (solid 

lines);  and  theTecond-order  wave  arising  from  rescatter 

of  the  first-order  wave  off  the  second  heterogeneity  (dotted  lines). 


p2  and  separation  between  the  two  objects  pi-p2  (see 
Fig.  2).  From  Eq.  (2),  expressions  can  be  written  for 
measurements  of  4>*i  and  4^2  made  at  detector  posi¬ 
tion  pd  in  the  presence  of  single  objects: 

4>*l(pd)  =  d>inc(Pd)  +  <J>acat>in’*1(Pd)» 

<t>*2(Pd)  =  ^.nc(Pd)  +  ^scaU2"_1(pd),  (4) 

where  kl  denotes  the  presence  of  the  first  body  alone 
and  k2  denotes  the  presence  of  the  second  body. 
With  simple  algebraic  manipulation  of  Eqs.  (2H4)> 


k=l  k-2 

Fig.  2.  Schematic  detailing  the  geometry  used  in  the  calculation 
of  scattered  waves  from  analytical  expression.  The  centroid  of  the 
cylinder  is  the  origin,  with  z  denoting  the  length,  angle  d  denoting 
the  angle  in  the  plane  containing  the  source  and  the  detector,  and 
r  denoting  the  radial  direction. 

predicted  by  measurements  of  d^ifpd).  d>«(P(f)>  an(1 

4>inc(Pd): 


d>*U2(pd)  =  ^*  =  l(P d)  +  <*\=2(Prf)  -  ^inc(pd)-  (6) 

Because  we  report  our  results  in  terms  of  phase  shift 
and  amplitude  demodulation,  Eq.  (6)  can  be  written 
as 


.,/tfn  sin  +  M* 2  sin  0*2  - 

^2(Pd)  tan  cos  +  Mk2  cos  9*2  - 


M.nc  sin  0i, 


(7) 

(8) 


an  expression  accounting  for  higher-order-scattering 
effects  can  be  written  for  the  photon  density  wave  in 
the  presence  of  both  objects  <^*ij*2(Pd)' 

4)*U2(Pd)  =  ^lfPd)  +  $*2(P</)  “  ^inc(pd) 

+  4>3Cat,*l,*2  (Pd)  4"  4>scat.*l>2  (Pd)  "*"•••• 

(5) 


B.  Analytical  Predictions  of  Phase  Shift  and  Amplitude 
Modulation  Owing  to  Two  Perfect  Absorbing  Cylinders 

The  complex  fluence  describing  the  propagation  of  a 
photon  density  wave  can  also  be  computed  analyti¬ 
cally  from  the  Helmholtz  equations.15-16  The  com¬ 
plex  fluence  propagating  from  a  source  point  at  ps 
through  an  infinite  random  medium  and  received  at 
position  p,  in  the  absence  of  an  optical  heterogeneity, 
is  given  by 


For  this  study,  it  is  assumed  that  third-  and  higher- 
order  scattered  waves  are  considered  to  have  insig¬ 
nificant  contributions  to  the  measured  fluence  when 
compared  with  those  of  first-  and  second-order  scat¬ 
tered  waves. 

From  frequency  domain  measurements  of  phase 
shift  0  and  ac  amplitude  M,  values  of  the  complex 
fluence  $  =  Afr-  ^  can  be  obtained  (1)  in  the  presence 
of  the  first  object  alone,  4>*i(prf);  (2)  in  the  presence  of 
the  second  object  alone,  4>*2(Pd)>  (3)  in  the  presence  of 
both  objects,  4>*i>2(Pd);  311(1  (4) in  the  absence  of  any 
inhomogeneities,  If  second-order  effects, 

i.e.,  4>SCat.*W2n=2(Pd).  are  negligible,  from  Eq.  (5)  it 
follows  that  measurements  of  d^i^Prf)  should  be 


exp 


d^incfp)  ^source(Ps) 


.(-  |X0  +  iw/Cn 


1/2 


(P  “  Pd) 


4irZ)(p  -  pd) 


(9) 


where  Scarce  describes  the  phase  and  strength  of 
modulation  of  the  source  located  at  position  ps;  c„  is 
the  speed  of  light  in  the  medium;  and  D  is  the  optical 
diffusion  coefficient  that  is  governed  by  the  absorp¬ 
tion  p.a  and  the  isotropic  scattering  |x's,  coefficients  of 
the  continuous  or  homogeneous  medium: 

D  =  i/3(m  + 1 V).  ‘  (10) 


9060  APPLIED  OPTICS  /  Vol.  36,  No.  34  /  1  December  1997 


In  this  study,  we  approximate  a  point-modulated 
source  at  the  periphery  of  a  large  cylinder  as  a  point- 
modulated  source  located  at  the  interface  of  a  semi¬ 
infinite  medium.  In  this  case  the  fluence  is  written 
in  cylindrical  coordinates  and  at  the  surface  (z  =  0) 
given  by17 


^inc(Prf)  ^aourceCPa) 


4'irD 


exp 


~D~ 


1/2 


[(ps  “  Pdf  +  (z  "  Zq )2] 


(1/2 


exp 


[(p,  -  Pd)2  +  (2  -  Zo)2] 
\  1/2 

M-at-  iiii/cT 


1/2 


D 


[(ps  -  Pd)2  +  (z  +  Zo)2] 


1/2 


[(Ps  ~  Pd)2  +  (Z  +  Zq)2] 


1/2 


(11) 


where  z0  is  one  isotropic  scattering  length  (l/p.'s). 

The  first-order  (n  =  1)  scattered  wave  from  the  kl 
infinitely  long  cylinder  in  an  infinite  medium  is  given 

byl5.18 


^acauU  (Pd)  ^inc 


cos(m-&)cos(pz) 


xKn 


2  .  |  -V-a+iu/Cn 

p  +  l — 3 - 

o  i  -p-o  ici )/cn\ 
P  + 


1 1/2 


Pd 


D 


1/2 


DxIm’(x)Im(y)  -  Dkl’yIJx)Im'(y) 


DxKm'(x)Im(y )  -  Dkl'yKJx)Im'(y) J 


dp, 


(12) 

where  Dk'  is  the  optical  diffusion  coefficient  within 
the  cylinder  k;  Im  and  Km  are  modified  Bessel  func¬ 
tions;  and  parameters  x  and  y  are  given  as 


x  = 


y  = 


-H„  +  tw/Cq 

D 


2  |  ~P-0' +  ito/c, 

P  + 


Dk' 


1/2 

a*. 

11/2 

a*; 


and  ak  is  the  radius  of  cylinder  k  1.  Radius  pd,  angle 
-d,  and  length  2  denote  the  coordinates  of  the  detector 
relative  to  the  center  of  the  infinite  cylinder  kl  (Fig. 
2).  The  incident  wave  upon  the  infinite  cylinder  4>inc 
is  computed  from  Eq.  ( 1 1).  A  similar  expression  can 
be  written  for  $acat,*2',“1. 

We  calculated  second-order-scattering  contribu¬ 
tions  by  evaluating  the  scattered  wave  originating 


from  the  second  object  {, k2 )  as  the  incident  wave  upon 
the  first  object  (£1).  In  other  words, 


^acaUl"  2(Prf)  —  *^»ait*2 


"if 

m  =  l 


cos(mfl)cos(pz) 


xK, 


P2  + 


-IXa  +  Iw/Cn 


\'2 


Pd 


xKn 


2  ,  /-  M-a  +  Wc, 
P  + 


D 


11/2 


'  DxIm'(x)Im(y)  -  Dkl'yIm(x)Im'(y) 
DxKm'(x)Im(y)  -  Dkl'yKm(x)Im\y)\ 


dp, 


(13) 


where  the  incident  wave  upon  cylinder  k  1  is  now  the 
first-order  scattered  wave  <b?Cat.A2n=1  that  is  com¬ 
puted  from  Eq.  (12).  A  similar  expression  can  be 
written  for  d>8Catjk2n=2-  Because  the  optical  proper¬ 
ties  and  the  diameter  of  the  two  cylinders  were  iden¬ 
tical  in  this  study,  we  consider  the  case  in  which  Dx' 
=  D2'  and  az  =  a2. 

With  the  approach  described  above,  the  fluence  de¬ 
tected  at  position  p<*  can  be  computed  inclusive  of 
first-order-scattering  effects  [i.e.,  <!>(pd)  =  d’inJ.pJ  + 
2*=i  <t)Scat/,=:i(Pd)]>  as  wel1  as  second-order  scatter¬ 
ing  effects  [i.e.,  d>(Prf)  =  d)inc(prf)  +  S„=1  -*=1 
d)scatjtn(Pd)].  From  final  values  of  complex  fluence, 
the  phase  shift  and  amplitude  demodulation  can  be 
predicted  from  the  simple  relationships: 


6(p d)  =  tan 


_t  Im[d>(Prf)] 
Re[«>(pd)] 


M{Pd)  =  ({Im[d>(pi)]}2  +  {Re[<t>(pd)]}2)1/2. 


(14) 

(15) 


3.  Materials  and  Methods 


A.  Experimental  Measurements  of  Phase  Shift  Owing  to 
Two  Cylindrical  Absorbers 

To  experimentally  measure  MklJt2(pd),  Mkl(pd), 

Mk2(p d),  Afinc(prf).  0*lJfe2(prf)>  9*l(Prf)>  0*2(Pd)> 

einc(pd),  frequency  domain  measurements  were  made 
with  an  apparatus  that  uses  picosecond  pulsed  light 
at  780  nm  with  an  average  power  of  1.3  W.  Details 
of  the  apparatus  are  described  elsewhere.19  The 
phantom  consisted  of  a  plexiglass  cylinder  (with  a 
16.5-cm  diameter  and  20-cm  height)  filled  with  a 
0.5%  Intralipid  solution  (Kabi  Pharmacia,  Inc.,  Clay¬ 
ton,  N.C.).  As  illustrated  in  Fig.  3,  light  was  deliv¬ 
ered  to  a  peripheral  point  on  the  cylinder  with  a 
1000-p.m  fiber  (Model  HCP-M1000T-08,  Spectron 
Specialty  Optics  Co.,  Avon,  Conn.)  and  collected 
through  a  second  fiber  located  4  circumferential  cm 
from  the  incident  source.  Bakelite  plastic  rods  (di¬ 
ameter,  3.175  mm)  were  painted  black  to  provide 
perfectly  light-absorbing  inhomogeneities.  Mea¬ 
surements  of  0  ( )  and  M(pd)  at  80  MHz  were  con¬ 
ducted  as  the  rods  were  moved  in  tandem  along  the 
plane  perpendicular  to  a  line  connecting  the  source 
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measurements. 


and  detector.  Positioning  was  achieved  with  a  mo¬ 
tion  controller  (Model  PMC200-P,  Newport  Corp.,  Ir¬ 
vine  Calif.)  and  a  motorized  actuator  (Newport 
Model  850B).  The  actuator  position  was  accurate  to 
within  0.0050  mm.  Phase  and  ac  modulation  were 
recorded  as  the  distance  between  the  center  of  the 
first  absorbing  cylinder,  and  the  wall  of  the  phantom 
was  varied  from  0  to  5  cm  in  20  increments.  When 
the  two  objects  were  moved  in  tandem,  their  separa¬ 
tion  distances  (Pl-p2)  between  the  centers  of  the  two 
perfect  absorbers  were  6,  10,  and  20  mm.  Three 
measurements  were  taken  at  each  position.  Mea¬ 
surements  of  phase  shift  and  ac  amplitude  demodu¬ 
lation  are  reported  relative  to  the  absence  case  or  9inc 
or  M;nr. 


B.  Analytical  Prediction  of  Phase  Shift  and  Amplitude 
Modulation  Owing  to  Two  Absorbing  Cylinders 
In  addition  to  the  measurements  of  second-order  in¬ 
teractions  with  the  experimental  approach  described 
above,  predictions  of  interactions  between  absorbers 
were  computed  analytically  with  Eqs.  (11M15)  with 
a  modified  version  of  an  algorithm  written  by  Boas 
that  is  available  through  the  internet.20  Complex 
fluence  incorporating  first-  and  second-order¬ 
scattering  contributions  was  used  to  calculate  phase 
shift  owing  to  two  infinite  cylinders  that  effectively 
act  as  perfect  absorbers  (p.a  =  2  cm  1,  m  =  10  cm  ) 
in  a  turbid,  semi-infinite  medium  mimicking  our 
phantom  (fiQ  =  0.02  cm"1,  (x/  =  10  cm"  ).  Zero 
fluence  boundary  conditions  were  assumed  and  used 
in  the  algorithm. 


C.  Numerical  Prediction  of  Phase  Shift  and  Amplitude 
Demodulation  Owing  to  Two  Absorbing  Heterogeneities 
In  addition  to  experimental  measurements  and  ana¬ 
lytical  calculations,  two-dimensional  finite-element 
computations  were  conducted  to  predict  the  phase 
shift  and  amplitude  demodulation  resulting  from 
first-  and  second-order  interactions.  The  computa¬ 
tions  were  performed  on  an  Ultra  2  Sun  Sparc  work¬ 


station  with  a  matlab  partial  differential  equation 
!£  We  solved  the  frequency  domain  diffusion 
equation  for  light  propagation  to  determine  the  flu- 
ence  from  an  80-MHz  sinusoidally  modulated  light 
source  The  simulated  phantom  was  an  infinite  cyl¬ 
inder  16.5  cm  in  diameter.  The  optical  properties  of 
the  medium  were  set  to  mimic  the  experLmental  con- 
ditions  of  0.02  cm  for  the  absorption  coefficient  ^ 
and  10  cm-1  for  the  isotropic  scattering  coefficient 
ll  '  We  modeled  the  nearly  perfect  absorbers  m 
these  computations  as  infinite  cylindrical  rods  and 
we  approximated  by  increasing  the  absorption  coef¬ 
ficient  to  100  times  that  of  the  surrounding  medium 
(„  =20  cm-1).  The  scattering  coefficient  tor^the 

object  was  set  to  that  of  the  surroundings  W  =  10 


w do  act  w  - -  “  .  '  —v 

cm-)  to  mimic  the  analytical  computations  The 


--1\ 


cm  )  WJ  limuiu  WAV,  -  *  ,  y 

phantom  was  discretized  into  66,048  triangular  ele¬ 
ments  that  contained  33,281  nodes.  A  partial  cur¬ 
rent  boundary  condition  was  used  to  approximate 
light  reflection  and  transmission  across  the  boundary 
of  the  phantom  as  would  be  expected  m  the  experi¬ 
mental  measurements.  Although  our  meshing  did 
not  permit  representation  of  the  perfect  absorber  as  a 
volume  excluded  for  light  transport,  we  approxi¬ 
mated  the  heterogeneity  with  a  high  absorption  co¬ 
efficient  because  others  have  shown  that  the 
propagation  characteristics  are  comparable.  It  is 
noteworthy  that,  although  these  computations  do  not 
exactly  mimic  the  experimental  measurements  de¬ 
scribed  below,  they  nonetheless  adequately  describe 
the  contributions  of  second-order  scattering  effects. 
The  forward  solution  was  obtained  for  each  absorbing 
cylinder  alone  and  then  in  combination.  The  results 
were  analyzed  in  a  manner  similar  to  the  analysis  of 
the  experimental  results  with  Eq.  (7)  and  (8). 


4.  Results  and  Discussion 


A.  Experimental  Results 

Figures  4(a)-4(c)  and  5(a)-5(c)  illustrate  the  phase 
shift  and  the  amplitude  demodulation  changes  as  a 
function  of  position  of  the  pair  of  perfect  absorbers 
whose  centers  are  separated  by  6,  10,  and  20  mm. 
The  symbols  denote  individual  measurements  of 
phase  shift  difference  in  the  presence  of  the  two  hid¬ 
den  objects,  and  the  symbols  connected  by  the  lines 
denote  the  phase  shift  and  the  amplitude  demodula¬ 
tion  calculated  from  Eq.  (7)  and  (8)  with  measure¬ 
ments  made  in  the  absence  and  in  the  presence  of  a 
single  individual  absorber.  The  measured  values  of 
phase  shift  and  amplitude  demodulation  denoted  by 
the  open  symbols  are  therefore  reflective  of  higher- 
order  contributions,  whereas  the  solid  line  reflects 
only  first-order  contributions.  Paired  Student’s  t- 
test  shows  that  there  is  a  significant  difference  (p  < 
0.005)  between  the  set  of  phase  shift  and  amplitude 
demodulation  measurements  and  that  predicted  by 
Eqs.  (7)  and  (8)  for  two  perfect  light-absorbing  objects 
separated  by  6,  10,  and  20  mm  and  positioned  at 
various  distances  from  the  source  and  the  detector 
[Figs.  4(a)-5(c)].  Our  results  also. show  that  the 
magnitude  of  second-order  effects  is  greatest  for  two 
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Fig.  4.  Experimental  values  of  6*1>2-^mc  the  phase  difference  (in 
degrees)  relative  to  the  absence  condition  as  a  function  of  distance 
from  the  source- detector  pair.  The  symbols  denote  individual 
measurements  in  the  presence  of  two  cylinders  separated  by  dis¬ 
tances  of  (a)  6,  (b)  10,  and  (c)  20  mm,  whereas  the  line  connects 
predictions  from  Eq.  (7)  and  measurements  of  0*1(  0*2,  and  0inc.  The 
error  bars  denote  the  propagation  of  measurement  errors  (standard 
deviation)  associated  with  0A1,  0*2.  and  0^.  The  x  axis  is  reported 
as  the  distance  between  the  wall  and  the  first  cylinder  (61).  The 
asterisks  denote  significant  difference  (p  <  0.005,  paired  Student’s 
t-test)  between  the  values  experimentally  measured  and  those  ob¬ 
tained  from  Eq.  (7). 


absorbing  heterogeneities  spaced  6  mm  apart  and 
becomes  smaller  at  10  and  20  mm.  The  paired  £-test 
indicates  that  there  are  significant  differences  at  the 
99.5%  confidence  levels  (indicated  by  the  asterisks  in 


Distance,  p,  from  source-detector  pair  [cm] 
(C) 


Fig.  5.  Experimental  values  of  MklJt2/Minc  (in  arbitrary  units) 
relative  to  the  absence  condition  as  a  function  of  distance  from  the 
source- detector  pair.  The  symbols  denote  individual  measure¬ 
ments  in  the  presence  of  two  cylinders  separated  by  distances  of  (a) 
6  (b)  10,  and  (c)  20  mm,  whereas  the  line  connects  predictions  from 
Eq.  (8)  and  measurements  of  Mklf  Mk2t  and  Afinc.  The  error  bars 
denote  the  propagation  of  measurement  errors  (standard  devia¬ 
tion)  associated  with  Mkl,  Mk2,  and  Minc.  The  x  axis  is  reported  as 
the  distance  between  the  wall  and  the  first  cylinder  (fcl). 


Fig.  4)  between  individual  experimental  phase-shift 
values  that  reflect  first-  and  higher-order  scattering 
contributions  to  the  detected  signal  and  those  com¬ 
puted  values  that  are  indicative  of  first-order  effects 
only.  From  the  data  for  6-mm  absorber  separation 
illustrated  in  Fig.  4(a),  it  can  be  seen  that  the  actual 
experimental  measurements  of  phase  shift  change 
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Fig.  6.  Values  of  6*1>2-0,„c  (in  degrees)  calculated  from  the  ana¬ 
lytical  prediction  of  first-order  (solid  line)  and  inclusive  of  second- 
order  (dashed  line)-scattering  effects  as  a  function  of  distance  (cmi 
from  the  source-detector  pair  for  two  absorbing  cylinders  sepa¬ 
rated  by  (a)  6,  (b)  10,  and  (c)  20  mm.  The  *  axis  is  reported  as  the 
distance  between  the  wall  and  the  first  cylinder  (&1)  and  the  phase 
shift  reported  relative  to  the  absence  case. 


caused  by  two  objects  are  smaller  than  those  pre¬ 
dicted  by  Eq.  (7),  in  which  second-  and  higher-order 
perturbations  are  not  accounted  for.  At  greater  dis¬ 
tances  from  the  source  and  detector,  agreement  be¬ 
tween  experimental  and  predicted  phase  shift  change 
suggests  that  second-order  effects  may  indeed  be  neg¬ 
ligible,  even  for  the  smallest  separation,  but  only  in  a 
region  where  the  objects’  contributions  to  the  de¬ 
tected  signal  is  comparatively  small.  This  is  reason¬ 


able  because  second-order  effects  are  expected  to 
increase  with  proximity  to  the  source  and  the  detec¬ 
tor.  From  our  studies,  it  appears  that  separations 
greater  than  20  mm  are  necessary  for  their  accurate 
resolution  through  perturbative  reconstruction  ap¬ 
proaches  when  perfect  absorbers  are  involved. 

B.  Analytical  Computations 

Our  experimental  results  are  also  validated  by  the 
analytical  predictions  that  account  for  first-  and 
second-order  scattering  contributions.  Figures 
6(a)-6(c)  depict  the  qualitative  trends  seen  m  the 
experimental  data  presented  for  phase  difference  in 
Figs  4(a)-4(c).  There  appears  to  be  good  agreement 
between  the  trends  observed  experimentally  with 
various  separation  distances  and  analytical  predic¬ 
tions.  Specifically,  higher-order  contributions  sig¬ 
nificantly  perturb  the  detected  phase  shift  values 
when  the  separations  between  the  center  of  two  ab¬ 
sorbing  cylinders  are  6  mm  [Fig.  6(a)]  and  10  mm 
[Fig.  6(b)].  The  contribution  of  higher-order  scatter¬ 
ing  from  cylindrical  absorbers  is  not  significant  when 
the  separation  distance  is  20  mm  [Fig.  6(c)].  How¬ 
ever  comparison  with  experimental  results  shows 
that  there  are  differences  in  the  absolute  magnitude 
and  shape  of  the  phase  shift  change  versus  object 
distance.  These  distances  are  most  likely  caused  by 
the  differences  in  the  boundary  conditions,  geometry, 
and  absorber  strength  between  the  analytical  and  the 
experimental  results.  Nonetheless,  the  trends  con¬ 
firm  experimental  results  that  the  presence  of 
second-order  perturbations  reduces  the  phase  shift 
change  owing  to  two  light-absorbing  bodies.  Ne¬ 
glecting  second-order  effects  could  cause  an  underes¬ 
timation  of  absorption  strength  or  size  when 
reconstruction  algorithms  based  on  first-order  per¬ 
turbations  are  used. 

In  addition,  we  investigated  the  variation  of 
second-order  perturbations  with  modulation  fre¬ 
quency  as  shown  in  Fig.  7.  The  simulated  data  of 
phase  shift  change  versus  position  of  the  heterogene¬ 
ities  is  depicted  for  two  cylindrical  absorbers  (d  = 
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Fig.  8.  Finite  element  computations  of  6*i>2*einC  (in  degrees) 
relative  to  an  absence  condition  for  considering  only  first-order 
(solid  line)  and  including  second-order  (dashed  line)  perturbations 
as  a  function  of  distance  (cm)  from  the  source- detector  pair  for 
absorbing  cylinders  separated  by  (a)  6,  (b)  10,  and  (c)  20  mm.  The 
x  axis  is  reported  as  the  distance  between  the  wall  and  the  first 
cylinder  (kl)  and  the  phase  shift  reported  relative  to  the  absence 
case. 


3.125  mm,  |xa  =  2  cm-1)  separated  by  6  mm  in  a 
semi-infinite  medium.  Again,  the  phase  shift  is  re¬ 
ported  relative  to  the  absence  case,  and  the  distance 
is  reported  from  the  wall  to  the  center  of  the  first 
cylinder.  Three  frequencies  were  evaluated:  80, 
160,  and  240  MHz.  The  absolute  value  of  phase 
shift  increases  with  frequency,  but  the  contribution  of 
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Fig.  9.  Finite  element  computations  of  MklJk2/Minc  (hi  arbitrary 
units)  referenced  to  an  absence  condition  for  both  the  first-order 
(solid  line)  and  the  second-order  (dashed  line)  perturbations  as  a 
function  of  distance  (cm)  from  the  source- detector  pair  for  absorb¬ 
ing  cylinders  separated  by  (a)  6,  (b)  10,  and  (c)  20  mm.  The  x  axis 
is  reported  as  the  distance  between  the  wall  and  the  first  cylinder 
(Al)  and  the  modulation  ratio  is  reported  relative  to  the  absence 
case. 


second-order  effect  remains  roughly  the  same  abso¬ 
lute  value.  Consequently  the  error  in  assuming  neg¬ 
ligible  second-order  scattering  effects  becomes 
smaller  at  increasing  frequencies.  This  is  expected 
because  increased  damping  of  a  rescattered  photon 
density  wave  occurs  at  higher  modulation  frequen¬ 
cies.  While  these  results  intimately  depend  upon 
the  choice  of  optical  properties,  they  nonetheless 
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point  out  that  the  error  in  negleding  second-order 
scattering  contributions  in  analytically  based  recon¬ 
structions  becomes  smaller  with  increasing  modula¬ 
tion  frequencies.  Of  course  the  error  reduction 
occurs  at  the  expense  of  interrogating  a  smaller  vol¬ 
ume  of  tissue  with  a  source  modulated  at  an  in¬ 
creased  frequency.21 

C.  Finite  Element  Computations 

Our  experimental  results  were  validated  not  only  by 
analytical  predictions  but  also  with  finite  element  com¬ 
putations.  Figures  8(a)-9(c)  show  the  numerical 
computations  that  correspond  to  the  experimental 
phase  shift  data  presented  in  Figs.  4(a)-5(c).  Again 
there  appears  to  be  good  agreement  between  the 
trends  of  the  experimental  and  the  analytical  results 
and  the  numerical  solutions.  Also,  higher-order  con¬ 
tributions  perturb  the  detected  phase  shift  values 
when  the  separation  distances  between  the  two  rods 
are  6  mm  [Fig.  8(a)]  and  10  mm  [Fig.  8(b)].  The  con¬ 
tribution  of  higher-order  scattering  from  the  cylindri¬ 
cal  absorbers  is  not  significant  when  the  separation 
distance  is  20  mm  [Fig.  8(c)].  Figures  9(a)-9(c)  show 
the  amplitude  demodulation  relative  to  an  absence 
condition.  The  modulation  data  also  show  results 
similar  to  those  obtained  experimentally.  These  com¬ 
putation  results  confirm  that  the  presence  of  second- 
order  perturbations  is  important  for  two  light- 
absorbing  objects  that  are  less  than  20  mm  apart. 

5.  Conclusions 

In  summary,  our  experimental  measurements  show 
that  the  contributions  of  higher-order  scattering  of 
propagating  photon  density  waves  may  not  always  be 
insignificant.  While  our  analytical  and  numerical 
computations  do  not  exactly  reproduce  experimental 
conditions  (i.e.,  two-dimensional  finite  element,  semi¬ 
infinite  geometry,  etc.),  they  nonetheless  demonstrate 
that  the  experimental  trends  can  be  attributed  to 
second-order  effects.  Under  conditions  of  high  con¬ 
trast  caused  by  absorption,  analytical  approaches  to 
the  inverse  imaging  algorithm  may  restrict  the  reso¬ 
lution  and  the  sensitivity  of  biomedical  optical  imaging 
performed  in  the  frequency  domain.  An  analogy  can 
also  be  drawn  for  time  domain  and  cw  reconstruction 
approaches  that  do  not  deploy  full  solutions  to  the 
diffusion  equation  to  account  for  the  interdependence 
of  voxel  optical  properties  on  measured  fluences.  Cer¬ 
tainly  our  results  are  based  on  the  worst-case  scenario 
of  perfect  absorbers  and  may  have  less  impact  on  the 
image  reconstructions  involving  imperfect  absorbers. 
Under  conditions  in  which  multiple  heterogeneities 
are  contrasted  with  their  surroundings  on  the  basis  of 
scattering  or  fluorescence,  the  first-order  perturbation 
assumption  may  not  be  as  restrictive.  Indeed,  if  the 
nonlinearity  associated  with  second-order  scattering 
effects  is  small,  noniterative  Bom  and  Rytov  approxi¬ 
mations  are  especially  attractive  because  image  recon¬ 
struction  is  not  computation  intensive.  Nonetheless, 
our  results  suggest  that  inverse  imaging  algorithms 
that  depend  solely  on  first-order  perturbations  caused 
by  local  changes  in  tissue  absorption  properties  may 


not  provide  a  reconstruction  as  accurate  as  those  algo¬ 
rithms  that  depend  on  the  full  numerical  calculation  of 
the  diffusion  equation  with  specified  boundary  condi¬ 
tions. 

6.  Appendix  A 

1.  Nomenclature 

a*  Radius  of  cylinder  k  (cm) 

cn  Speed  of  light  through  the  medium  (cm/s) 

D  Optical  diffusion  coefficient  of  homogeneous 

medium  (cm) 

Dk>  Diffusion  coefficient  inside  the  cylinder  (cm) 

In  Modified  Bessel  function 

In.  Derivative  of  modified  Bessel  function  In 
Kn  Modified  Bessel  function 

Kn.  Derivative  of  modified  Bessel  function  Kn 
M  ac  demodulation  of  the  incident  light 
(mW/cm2) 

m  Number  of  objects 

p  Integration  variable  in  Helmholtz  equation  de¬ 

scribing  scatter  from  an  infinite  cylinder 
(cm“l) 

Source  Strength  of  modulated  source  at  position  p„ 
represented  as  a  complex  number  of  amplitude 
and  phase  (mW) 

x  Variable  in  Helmholtz  equation  describing 

scatter  from  an  infinite  cylinder  (cm"1) 
y  Variable  in  Helmholtz  equation  describing 

scatter  from  an  infinite  cylinder  (cm 
z  Axial  direction  of  the  cylindrical  heterogene¬ 

ities  or  length  (cm) 

Z  Greek 

<P  Photon  fluence  (mW/cm2) 

<Pinc  Photon  fluence  of  incident  wave  or  in  the  ab¬ 
sence  of  heterogeneities  (mW/cm2) 

< Photon  fluence  arising  from  scattered  wave 
(mW/cm2) 

p  Position  vector  (cm) 

ps  Position  of  the  source  (cm) 

pd  Position  of  detector  (cm) 

p*  Center  position  of  object  k  (cm) 

0  Phase  shift  of  light  wave  (deg  or  rad) 

Angle  between  ps  and  pd 

\iQ  Absorption  coefficient  of  homogeneous  sur¬ 
roundings  (cm"1) 

V  Isotropic  scattering  of  homogeneous  surround¬ 
ings  (cm"1) 

V  Absorption  coefficient  of  cylinder  (cm" l) 

IV  Isotropic  scattering  of  cylinder  (cm"1) 

3.  Subscripts  and  Superscripts 

n  Order  of  perturbation  or  nth-order  scattering  ef¬ 
fect 

k  Index  denoting  cylindrical  object  one,  kl;  two,  k2 ; 
or  both  objects,  kl,  k2 

inc  Index  denoting  absence  measurement  or  condi¬ 
tion 
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Abstract:  The  development  of  non-invasive,  biomedical  optical  imaging  from  time- 
dependent  measurements  of  near-infrared  (NIR)  light  propagation  in  tissues  depends 
upon  two  crucial  advances:  (i)  the  instrumental  tools  to  enable  photon  “time-of- 
flight”  measurement  within  rapid  and  clinically  realistic  times,  and  (ii)  the 
computational  tools  enabling  the  reconstruction  of  interior  tissue  optical  property 
maps  from  exterior  measurements  of  photon  “time-of-flight”  or  photon  migration.  In 
this  contribution,  the  image  reconstruction  algorithm  is  formulated  as  an 
optimization  problem  in  which  an  interior  map  of  tissue  optical  properties  of 
absorption  and  fluorescence  lifetime  is  reconstructed  from  synthetically  generated 
exterior  measurements  of  frequency-domain  photon  migration  (FDPM).  The  inverse 
solution  is  accomplished  using  a  truncated.  Newton’s  method  with  trust  region  to 
match  synthetic  fluorescence  FDPM  measurements  with  that  predicted  by  the  finite 
element  prediction.  The  computational  overhead  and  error  associated  with 
computing  the  gradient  numerically  is  minimized  upon  using  modified  techniques  of 
reverse  automatic  differentiation. 

©1999  Optical  Society  of  America 

Keywords:  (110:6960)  biomedical  optical  imaging;  (170:3010)  image  reconstruction;  (9:999)  Galerkin 
finite  element  method;  (9:9999)  Truncated  Newton’s  method;  (9:9999)  reverse  automatic  differentiation, 
(9:9999)  optimization;  (100:3190)  inverse  problem;  (260:2510)  fluorescence 


References: 

1  D.  A.  Boas,  M.  A.  O’Leary,  B.  Chance  and  A.G.  Yodh,  “Scattering  of  diffuse  photon  density  waves  by  spherical 
heterogeneities  within  turbid  media:  analytic  solutions  and  applications,”  Proc.  Natl.  Acad.  Sci.,  91.  pp  4887-91, 
(1994). 

2  M.  A.  O’Leary.  D.  A.  Boas,  B.  Chance  and  A.G.  Yodh,  “Experimental  images  of  heterogeneous  turbid  media  by 
frequency-domain  diffusion  photon  tomogrpahy,”  Opt.  Lett.  20,  426-428,  (1995). 

3  R.  L.  Barbour.  H.  Graber,  Y.  Wang  J.  Chang  and  R.  Aronson,  “Perturbation  approach  for  optical  diffusion 
tomography  using  continuous-wave  and  time-resolved  data,”  in  Medical  Optical  Tomography:  Functional  Imaging 
and  Monitoring,  G.  Muller,  B.  Chance,  R.  Alfano,  J.  Beuthan,  E.  Gratton,  M.  Kashke,  B.  Masters,  S.  Svanberg  and  P. 
van  der  Zee,  eds.  (SPIE  Press,  Bellingham,  WA.,  1993),  pp  87-120. 

4  Y.  Yao,  Y.  Wang,  Y.  Pei,  W.  Zhu  and  R.L.  Barbour,  “Frequency-domain  optical  imaging  of  absorption  and 
scattering  by  a  Bom  iterative  method,"  J.  Opt.  Soc.  Am.  A,  14,  325-342,  (1997). 

5  W.  C.  Chew  and  Y.  M.  Wang,  “Reconstruction  of  two-dimensional  permittivity  distribution  using  the  distorted 
Bom  iterative  method,”  IEEE  Trans.  On  Medical  Imaging.  9,  218-225,  (1995). 

6  K.  D.  Paulsen  and  H.  Jiang,  “Spatially  varying  optical  property  reconstruction  using  a  finite  element  diffusion 
equation  approximation.”  Med.  Phys.  22,  691-701,  (1995). 

7  H.  Jiang,  K.  D.  Paulsen,  U.  L.  bsterberg,  B.W.  Pogue  and  M.  S.  Patterson.  “Optical  image  reconstruction  using 
frequency-domain  data  simulations  and  experiments,”  J.  Opt.  Soc.  Am.  A.  13,  253-266,  (1996). 

8  D.Y.  Paithankar.  A.  U.  Chen,  B.W.  Pogue,  M.  S.  Patterson  and  E.  M.  Sevick-Muraca.  "Imaging  of  fluorescent 
yield  and  lifetime  from  multiply  scattered  light  re-emitted  from  tissues  and  other  random  media,"  Appl.  Opt ,  36, 


2260-2272,  (1997). 

9  M.  Schweiger,  S.  R.  Arridge.  D.  T.  Delpy,  "  Application  of  the  finite-element  method  for  the  forward  and  inverse 
models  in  optical  tomography,"  J.  Math.  Img.  and  Vision  3,  263-283,  (1993) 

10  J.  J.  McKeown  ,"On  algorithms  for  sums  of  squares  problems,"  Towards  global  optimization,  edited  by  Dixon.  L. 
C.  W.  and  Szeeg,  G.  P.,  North-Holland  Amsterdam,  Holland,  (1975). 

1 1  M.  T.  Vespucci,  "An  efficient  code  for  the  minimization  of  highly  nonlinear  and  large  residual  least  squares 
functions,"  Optimization  18,  825-855,  (1987). 

12  R.  R.  Meyer,  "Theoretical  and  computational  aspects  of  nonlinear  regression,"  Nonlinear  Programming,  edited 
by  Rosen,  J.  B.,  Mangasarian,  O.  L.  and  Ritter.  K.,  Academic  Press,  New  York,  (1970). 

13  L.  B.  Rail,  Automatic  differentiation:  Techniques  and  application.  Lecture  notes  in  computer  science.  120, 
Springer  Verlag,  (1981). 

14  A.  Griewank,  "On  automatic  differentiation."  edited  Iri,  M.  and  Tanaka,  K.,  Mathematical  programming:  Recent 
developments  and  application,  Kluwer  Academic  Publishers,  pp  83-108.  (1989). 

15  T.  L.  Troy,  D.  L.  Page  and  E.  M.  Sev ick-Muraca,  "Optical  properties  of  normal  and  diseased  breast  tissues: 
prognosis  for  optical  mammography,"  J.  Biomed  Opt.  1,  342-355,  (1996). 

16  E.  M.  Sevick-Muraca,  G.  Lopez.,  T.  L.  Troy,  J.  S.  Reynolds  and  C.  L.  Hutchinson,  "Fluorescence  and  absorption 
contrast  mechanisms  for  biomedical  optical  imaging  using  frequency-domain  techniques,"  Photochemi  and 
Photobio l,  66:  55-64,  (1997). 

17  M.  A.  O’Leary,  D.  A.  Boas,  B.  Chance  and  A.G.  Yodh,  “Fluorescence  lifetime  imaging  in  turbid  media,”  Opt. 
Lett.. 21,  158-160,(1996). 

18  J.  Chang,  H.  L.  Graber  and  R.L.  Barbour,  “Luminescence  optical  tomography  of  dense  scattering  media.”  J.  Opt. 
Soc.  Am.  A.  14,  288-299,  (1997). 

19  T.  L.  Troy  and  E.  M.  Sevick-Muraca,  "Fluorescence  lifetime  imaging  and  spectroscopy  in  random  media,"  in 
Applied  Fluorescence  in  Chemistry \  Biology,  and  Medicine,  Rettig,  Strehmel,  Shrader,  Seifert,  eds.,  Springer  Verlag, 
pp.  3-36,  (1999). 

20  H.  Jiang,  “Frequency-domain  fluorescent  diffusion  tomography:  a  finite-element  based  algorithm  and 
simulations,”  Appl.  Opt.  37,  5337-5343,(1998). 

21.  Chang,  H.  L.  Graber  and  R.L.  Barbour.  “Improved  reconstruction  algorithm  for  luminescence  optical  tomography 
when  background  Iuminophore  is  present,”  Appl  Opt.  37 , 3547-3552.  (1998). 

22  J.  Lee  and  E.  M.  Sevick-Muraca,  “Lifetime  and  absorption  imaging  with  fluorescence  FDPM,”  Time-resolved 
fluorescence  spectroscopy  and  imaging  in  tissues,  E.  M.  Sevick-Muraca  (ed.).,  Proc.  Soc.  Photo-Opt.  Instrum.  Eng., 
3600:  (to  be  published),  (1999). 

23  A.  Ishimaru,  Wave  propagation  and  scattering  in  random  media,(  Academic  Press,  New  York,  1978). 

24  M.  Schweiger.  S.  R.  Arridge,  M.  Hiraka  D.  T.  Delpy,  "  The  finite-element  method  for  the  propagation  of  light  in 
scattering  media-  boundary  and  source  conditions,"  Med.  PJtys  22,  1779-1792,  (1995) 

25  R.  A.  J.  Groenhuis,  H.  A.  Ferwerda  and  J.  J.  Ten  Bosch,  "Scattering  and  absorption  of  turbid  material  determined 
from  reflection  measurements,"  Appl.  Opt.  22,  2456-2462,  (1983). 

26  O.  C.  Zienkiewcz  and  R.  L.  Taylor , The  finite  element  methods  in  engineering  science.  (McGraw-Hill.  New  York, 
1989). 

27  L.  C.  W.  Dixon  and  R.  C.  Price,  "Numerical  experience  with  the  truncated  Newton  method  for  unconstrained 
optimization,"  JOTA,  56,  245-255,  (1988). 

28  R.  Roy,  Image  reconstruction  from  light  measurements  on  biological  tissue,  Ph.  D.  thesis,  University  of 
Hertfordshire,  England, (1996). 

29  R.  S.  Dembo  and  T.  Steihaug,  "Truncated  Newton  algorithms  for  large-scale  unconstrained  optimization,"  Math 
Programming  26,190-212,(1983). 

30  R.  C.  Price,  Sparse  matrix  optimization  using  automatic  differentiation,  Ph.  D.  thesis.  University  of  Hertfordshire. 
U.  K.,  (1987). 

31  L.  Armijo  "Minimization  of  functions  having  Lipschitz  continuous  first  partial  derivatives."  Pacific  J. 
Mathematics  16,  1-3.  (1966). 

32  P.  Wolfe,  "Convergence  condition  for  ascent  method."  SIAM  Rev.,  11  226-253,  (1969). 

33  B.,  Christianson,  A.  J.,  Davies,  L.  C.  W.  Dixon,  R.  Roy  and  P.  van  der  Zee,  "Giving  reverse  differentiation  a 
helping  hand,"  Opt.  Meth.  And  Software  8,  53-67,  (1997). 

34  A.  J.,  Davies,  B.  Christianson,,  L.  C.  W.  Dixon,,  R.  Roy  and  P.  van  der  Zee,  "Reverse  differentiation  and  the 
inverse  diffusion  problem,"  Adv.  In  Eng.  Software  28,  217-221,  (1997). 

35  R.  E.  Wengert.  "A  simple  automatic  derivative  evaluation  program,"  Comm.  A.  C.  M.,  7,  463-464,  (1964). 


1,  Introduction 

Conventional  imaging  modalities  such  as  magnetic  resonance  imaging  (MRI)  and  x-ray 
computer-aided  tomography  (x-ray  CT)  provide  high-resolution  medical  imaging  enabled  by 
the  direct  geometric  correlation  between  incident  and  detected  radiation.  Yet  the  high  cost  of 
operating  hospital  MRI  facilities,  and  the  inability  to  detect  important  diseases  through  x-ray 


CT  imaging  suggests  an  opportunity  for  the  development  of  near-infrared  (NIR)  biomedical 
optical  imaging.  NIR  biomedical  optical  imaging,  or  optical  tomography,  depends  upon  the 
low  absorbance,  yet  high  scattering  of  non-mutagenic  near-infrared  light  in  the  “therapeutic 
wavelength  window”  (600-1000  nm)  enabling  it  to  safely  propagate  through  several 
centimeters  of  tissues.  While  the  propagation  of  low-energy  NIR  light  can  occur  safely,  the 
high  scattering  properties  of  tissues  renders  the  direct  geometric  correlation  between  incident 
and  detected  irradiation  invalid.  The  forward  imaging  problem,  (i.e.,  prediction  of  the 
propagation  of  NIR  light  through  tissues  when  a  map  of  the  tissue  optical  properties  is 
known)  can  be  described  in  terms  of  the  diffusion  approximation  to  radiative  transfer.  The 
inverse  imaging  problem  (i.e.,  prediction  of  the  interior  optical  properties  from  measurements 
of  light  propagation  made  at  the  exterior  tissue-air  interface  positions)  requires  solution  of  a 
series  of  equations  that  are  non-linear  in  the  optical  properties  to  be  estimated.  In  this 
contribution,  we  describe  a  novel  algorithm  to  estimate  a  solution  to  the  nonlinear,  inverse 
imaging  problem  for  fluorescence  frequency-domain  photon  migration  (FDPM).  In  the 
following  sections  of  this  introduction,  we  briefly  (i)  present  the  background  of  frequency- 
domain  photon  migration  (FDPM)  imaging,  (ii)  compare  prior  work  and  our  current  approach 
towards  solution  of  its  inverse  imaging  problem,  and  (iii)  introduce  the  concept  of  fluorescent 
contrast  agents  which  can  accelerate  the  convergence  of  the  inverse  imaging  solution. 

1.1  Frequency-domain  photon  migration:  forward  and  inverse  problems 

Frequency-domain  photon  migration,  depends  upon  launching  intensity  modulated  (30-200 
MHz)  light  at  the  interface  of  a  highly  scattering  medium  (such  as  tissue),  and  detecting  the 
intensity-modulated  wave  that  successfully  propagates  to  the  detector  located  a  distance  p 
away  from  the  incident  source.  Depending  upon  the  spatial  distribution  of  interior  absorption 
and  scattering  optical  properties,  the  detected  light  is  both  phase-delayed  by  0,  and  amplitude 
attenuated  by  factor  M  when  compared  to  the  incident  light  (see  Figure  1).  The  solution  to  the 
forward  imaging  problem  involves  predicting  the  phase-delay,  0,  and  amplitude 
demodulation,  M  as  a  function  of  angular  modulation  frequency,  co,  and  position,  r,  along  the 
tissue  surface  given  a  known  spatial  distribution  of  absorption  and  scattering  optical 
properties.  The  absorption  coefficient,  |Lia  ,  representing  the  inverse  mean  absorption  path 

and  the  isotropic  scattering  coefficient,  u.^ ,  equivalent  to  the  inverse  isotopic  scattering 

length,  govern  the  propagation  of  light  through  scattering  media,  such  as  tissues.  Since  near- 
infrared  light  is  multiply  scattered  in  tissues,  the  numerical  solution  to  the  diffusion 
approximation  of  the  radiative  transfer  equation  provides  solution  to  the  forward  biomedical 
optical  imaging  problem.  For  a  propagating  intensity  modulated  wave  of  light,  the  optical 
diffusion  equation  is  written  as 


—  V  •  (r )  VO  (r ,  co)]  + 


d>^.(r,co)  =0  on  SI  (l) 


where  represents  the  complex  number  describing  the  scalar  flux  of  photons  at  position  r; 
co  is  the  modulation  frequency;  c  is  the  speed  of  light  within  the  medium;  is  the 

absorption  due  to  chromophores  and  \ia  is  the  absorption  due  to  fluorophores 
(|i,^  =[ia  .  +l±aKf  )•  an  optical  diffusion  coefficient  equivalent  to  y|3(fj.^  +  M^v )j- 

From  numerical  solution  of  Eqn  (1)  with  a  Robin  boundary  condition,  the  complex  scalar  flux 
at  the  surface,  0^(=Af^,e),  and  the  phase  delay,  0,  and  amplitude  attenuation,  M,  can  be 


determined.  It  is  important  to  note  that  the  diffusion  equation  applies  when  ^  «  |l'!r  and 
Dx  is  approximated  by  )J. 


Figure  1  Schematic  of  fluorescence  photon  migration. 


The  successful  implementation  of  biomedical  optical  imaging  involves  the  effective 
solution  of  the  inverse  imaging  problem,  i.e.,  determining  the  interior  optical  property  map  of 
absorption  and  scattering  given  measurements  of  0  and  M  along  the  tissue  surface.  To  date, 
approaches  have  focused  upon  the  linearization  of  the  problem  using  first  order  Born  or 
Rytov  approximations  [1,  2].  Non-linear  optimization  employing  iterative  Born  [3,  4],  or 
distorted  Born  methods  [3,5]  require  accurate  information  regarding  the  normal  ’‘background” 
optical  properties  which  the  tissue  “heterogeneity”  resides.  Since  the  normal  tissue 
background  optical  properties  can  be  expected  to  vary  greatly,  this  a  priori  information  is  not 
realistically  feasible  for  biomedical  optical  imaging. 

Other  investigators  have  used  the  modified  Newton-Raphson  method  with  the 
Levenberg-Marquardt  regularization  as  the  central  inversion  step  [6,  7,  8  and  9].  While  not 
requiring  accurate  “background”  optical  properties  for  initial  parameter  estimate  guesses  for 
successful  reconstruction,  the  Newton-Raphson  with  Levenberg-Marquardt  regularization  has 
three  distinct  disadvantages  for  biomedical  optical  imaging.  The  first  disadvantage  is  that  the 
full  Hessian  matrix  is  not  fully  considered  making  the  method  inefficient  for  highly  non¬ 
linear  problems  such  as  optical  imaging.  Numerical  experiments  by  McKeown  [10]  and 
Veapucci  [111]  have  shown  that  while  the  Gauss-Newton  and  Levenberg-Marquardt  methods 
solve  nonlinear  least  squares  problems  efficiently  when  the  residuals  are  small  or  zero  at  the 
solution,  these  methods  are  less  efficient  than  the  quasi-Newton  method  when  the  residuals 
are  significantly  nonlinear.  Meyer  [12]  provides  theoretical  support  of  this  observation.  In  the 
work  presented  herein,  we  demonstrate  the  use  of  a  truncated  Newton  method  with  trust 
region  for  the  efficient  solution  of  the  large  scale,  non-linear  optical  imaging  problem.  The 
second  disadvantage  to  Newton-Raphson  approaches  is  that  the  methods  require  expensive 
numerical  calculation  of  the  Jacobian  that  can  be  subject  to  error.  In  this  work,  we  adapt 
modifications  of  the  reverse  automatic  differentiation  to  compute  the  gradient  defining  the 
search  direction  efficiently  and  without  numerical  error  [13,  14].  Finally,  the  Levenberg- 
Marquardt  method  requires  storage  of  the  full  Jacobian  matrix,  which  represents  a  limitation 
for  optical  imaging  in  three  dimensions  and  when  the  number  of  unknowns  increases.  Since 
the  truncated  Newton  method  requires  only  the  product  of  the  Hessian  with  the  direction  at 
any  point,  the  method  is  comparatively  more  efficient  for  large  problems.  Our  inversion 
strategy  is  presented  in  its  entirety  in  Section  3. 

While  the  difficulties  associated  with  the  nonlinear  and  large  scale  of  inverse 
imaging  problem  can  be  addressed  with  our  inversion  strategy,  additional  significant 
challenges  exist  when  attempting  to  reconstruct  absorption  and  scattering  properties  that  are 


physiologically  feasible.  In  a  study  involving  actual  optical  property  measurements  of  normal 
and  diseased  tissues,  Troy  et  al  [15]  found  that  the  optical  properties  of  absorption  and 
isotropic  scattering  of  normal  and  diseased  breast  tissues,  for  example,  may  not  provide 
sufficient  contrast  for  successful  image  reconstruction.  Indeed,  the  small  influence  of 
“heterogeneity”  absorption  and  scattering  on  FDPM  measurements  measured  at  the  boundary 
can  often  be  within  or  close  to  the  measurement  error  of  1  degree  in  phase  and  1%  in 
amplitude  attenuation. 

1.2  Fluorescence  frequency -domain  photon  migration:  forward  and  inverse  problems 

In  order  to  over  come  the  difficulties  associated  with  insufficient  contrast  and  costly  non¬ 
linear  optimization,  we  and  other  investigators  have  developed  fluorescence  FDPM  imaging 
approaches  and  inversion  strategies  [8,  16,  17,  18].  Fluorescence  FDPM  physically  depends 
upon  the  administration  of  a  fluorescent  contrast  agent,  or  another  agent  or  gene  vector,  which 
results  in  the  expression  of  a  fluorescent  protein  that  emits  in  the  near-infrared  wavelength 
regime.  When  activated  by  the  absorption  of  an  intensity  modulated  wave,  4>  v(r,  co),  an 
intensity  modulated  fluorescent  wave,  Om(r,co),  is  generated  at  the  same  modulation 
frequency,  but  is  amplitude  attenuated  and  phase-delayed  relative  to  the  activating  excitation 
wave.  The  phase-delay,  0m,  and  amplitude  attenuation,  Mm,  of  the  generated  fluorescent  wave 
can  be  as  high  as  90  degrees  and  a  factor  of  10’s  -100’s  of  the  incident  excitation  wave  owing 
to  the  fluorescence  properties  (lifetime  and  quantum  efficiency)  of  the  dye.  Specifically,  as 
the  biochemical  environment  of  the  fluorophore  changes,  its  lifetime,  x  (i.e.,  the  lifetime  of 
the  activated  fluorophore,  or  the  mean  time  between  absorption  of  excitation  and  release  ot 
fluorescent  light)  also  alters,  providing  discrimination  of  diseased  tissues  possibly  on  the 
basis  of  biochemical  environment.  The  detected  intensity  modulated  fluorescence  wave  that 
has  been  generated  within  and  which  has  propagated  to  the  tissue  boundary  possesses  a 
phase-delay,  0m,  and  amplitude  attenuation,  Mm,  that  is  exquisitely  more  sensitive  to 
embedded  “heterogeneities”  than  possible  with  absorption  FDPM  imaging  described  above 
[16]. 

Predictions  of  fluorescence  FDPM  measurements  of  phase-delay,  0m,  and  amplitude 
attenuation,  Mm,  are  achieved  through  the  solution  of  the  complex  fluorescence  fluence, 
Om(r,co),  at  the  boundary  from  the  fluorescence  optical  diffusion  equation: 


-V-[Dm(r)V<&m(r,co)]  + 


UX)  /_\ 

c  m 


Om(?,a))=  <()^a 


1 


1  -  ('cot 


<J>v(r,  co)  on  Q 


(2) 


where  Dm  is  the  optical  diffusion  equation  at  the  emission  wavelength 
(=  l/3|p.am  +^mD=l/3||^mJ;  liam  and  M-sm  are  the  absorption  and  isotropic  scattering 

coefficients  at  the  fluorescent  wavelength;  and  the  right  hand  term  describes  the  generation  of 
fluorescence  within  the  medium.  The  term  (p  represents  the  quantum  efficiency  of  the 
fluorescence  process  and  the  absorption  owing  to  fluorophore  is  represented  by  the  coefficient 
lu  .  Here  we  consider  first  order  relaxation  decays  only.  Note  that  the  source  term 

requires  coupling  with  the  solution  of  excitation  fluence  described  by  Eqn  (1),  where  by  the 
absorption  coefficient  at  the  excitation  light,  \ia  is  now  provided  by  the  naturally  occurring, 

non-fluorescent  chromophores,  \ia  .  and  the  fluorescent  contrast  agent,  .  The 

numerical  solution  for  the  excitation  fluence  distribution  Eqn  (1)  with  the  Robin  boundary 


condition  enables  prediction  of  from  Eqn  (2)  at  the  medium  boundary  and  determination 
of  0m  and  Mm.  It  is  important  to  note  that  the  diffusion  equation  applies  when  « JLl^  . 

With  the  added  unknowns  at  the  emission  wavelength,  the  disadvantages  using 
Newton-Raphson  optimization  for  fluorescence  FDPM  imaging  may  become  more  severe. 
Paithankar  et  al.,  [8]  and  Troy  and  Sevick-Muraca  [19]  employed  Newton-Raphson,  multi¬ 
grid  finite  difference  reconstructions  of  \*-ax_>m  an^  »  given  correct  priors  on  other 

optical  properties.  Recently,  Jiang  [20]  performed  similar  work,  employing  dual  meshing 
finite  element  methods  instead  of  finite  difference  methods.  Both  Paithankar  et  al.  [8]  and 
Jiang  [20]  resorted  to  removal  of  “spurious”  optical  property  values  or  “filtering”  using 
empirically  chosen  parameters  to  achieve  the  “correct”  image  distorted  through  inefficient 
numerical  computation  of  the  Jacobian.  O’Leary  et  al.,  [17]  and  Chang,  et  al,  [18]  employed 
Born  and  iterative  Born  which  again  not  only  required  known  “background”  optical 
properties,  but  which  also  was  unable  to  handle  fluorescence  in  the  “background.”  Later 
Chang  et  al,  [21]  demonstrated  the  ability  to  reconstruct  fluorescence  lifetime  under 
conditions  of  “imperfect”  uptake  given  the  spatial  distribution  of  absorption  owing  to 
fluorophore;  and  Lee  and  Sevick-Muraca  [22]  employed  an  iterative  Born-type  solution  of 
simultaneous  absorption  and  lifetime  from  fluorescence  FDPM  synthetic  data. 

1.3  Fluorescence  FDPM  imaging:  current  work 

In  this  contribution,  we  also  focus  upon  fluorescence  FDPM  imaging,  owing  to  the  significant 
influence  of  interior  fluorescent  heterogeneities  on  measured  signals  at  the  medium  boundary 
and  adapt  our  inversion  strategy  for  the  non-linear,  large  scale  optimization  required  for 
biomedical  optical  imaging  of  absorption  and  fluorescence  lifetime.  In  the  following  section, 
the  numerical  methodology  for  solving  the  FDPM  forward  imaging  problem  and  the 
simulator  used  for  generation  of  synthetic  data  sets  for  testing  our  inversion  algorithm  are 
presented.  Section  3  describes  our  formulation  of  the  inverse  imaging  problem  as  an 
optimization  problem  using  truncated  Newton’s  methods  with  trust  region  and  the  reverse 
automatic  differentiation  in  order  to  reduce  the  computational  overhead  associated  with 
optimization.  The  performance  of  our  algorithm  using  synthetic  data  sets  is  presented  in 
Section  4  and  the  conclusions,  prognosis  for  FDPM  imaging  and  ongoing  work  are 
summarized  in  Section  5. 

2.0  Model  prediction  of  FDPM 

2. 1  Boundary  conditions  for  FDPM  forward  simulator 

Prediction  of  FDPM  measurements  on  surface  dQ  resulting  from  diffusion  of  intensity- 
modulated  waves  within  volume  Q.  is  accomplished  from  solution  of  Eqn  (1)  and  (2)  using 
the  Robin  boundary  condition.  We  employ  the  Robin  boundary  condition  of  the  form  given 
by  Ishimaru[23,  24]: 

d>  v  (r,  co)—  2  y Dx  (r)  ^ — —  +  58  ( r  -  r5 )  =  0  on  dQ.  (3) 

3  n 

Where  5  is  the  complex  flux  density  associated  with  a  point  source  of  excitation  light 
located  at  rs  on  dQ;  n  is  the  unit  vector  normal  to  the  surface,  and  y  accounts  for  refractive 
index  mismatch  at  the  boundary  which  through  Snell’s  reflection,  results  in  internal  reflection 
at  the  surface. 

Groenhuis  et  al.  [25]  have  developed  an  empirical  approach  to  determine  the  value 
of  the  parameter  y  ,  viz. 


where 


y=(1  +  rj)/(1"r^) 


(4) 


r d  =  —  \AAnrei  +  0.72nr<?/  +  0.668  +  0.063//^/  (5) 

nrei  is  the  relative  index  of  refraction  of  the  scattering  medium  with  respect  to  the 
surrounding  medium. 

2.2  Galerkin  finite  element  formulation  for  solution  of<Px 

We  employ  the  Robin  boundary  condition  in  our  numerical  solution  of  Ox.m  employing  the 
Galerkin  finite  element  model  with  linear  triangular  elements.  In  our  formulation,  we  assume 
\ia  «[is  •  so  that  Dxm~\j3\iSxm  .  The  formulation  of  the  Galerkin  finite  element 

solution  of  <bx  begins  by  multiplying  Eqn  (1)  by  a  weighting  function,  Wj,  and  integrating  over 
the  domain  of  interest: 
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where  wj  {vv;- :  j  =  is  a  set  of  linearly  independent  weighting  functions.  N  is  the 

number  of  unknows. 

Using  Green's  theorem,  equation  (6)  becomes 
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(7) 


where  the  second  term  is  a  line  integral  along  the  external  boundary  and  Wj  is  assumed  to  be 
continuous  over  the  whole  domain.  Now  introducing  equation  (3)  in  the  line  integral 


=  -K<&x+s)wjdr  =  -\<&xWjdr  +  -jswjdr  (8) 

p  dn  Y  r  Y  p  Y  p 

Upon  combing  Eqns  (8)  and  (7),  once  can  write 
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Yr  Yr 


Suppose  that  the  solution  domain  Q,  is  divided  into  M  x  triangular  elements.  Now  we  can 
obtain  the  element  equation  from  Eqn  (9)  as  follows: 


(10) 
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Assuming  that  <t>eJ  ,DeJ and  (1^  vary  linearly  within  each  triangular  element  el  and  is 
constant,  we  may  then  write  <t>eJ ,DeJ and  [ieJ  for  each  element  as: 

axf 


<  =  iLj(*x)j  <ii> 

7=1 

where  the  L are  the  natural  co-ordinates  for  the  triangle  [26].  According  to  the  Bubnov- 
Galerkin  method,  the  weighting  functions  are  chosen  to  be  the  same  as  the  approximation 
functions  used  to  represent  <t>eJ  ,  that  is, 

Wj  =  Lj  for  j  =  1,2,..., N 


The  equation  (10)  can  be  written  as: 
M[  [  J  80  x1  3Lj  80  9Lj 
ei=i  nc  x  8x  8x  3y  3y 
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(12) 

After  integration  of  Eqn  (12)  and  rearrangement,  the  resulting  element  stiffness  equations 
become: 

|[Kf  +Ke^ +Ke3l]<t>e‘  =|re/  (13) 

el= 1  “  el 


Assembly  of  these  element  stiffness  matrices  K^and  local  vector  rel  yields  the  system 
equations  KOA.  =b  [26]. 

Since  5  denotes  the  complex  point  source  applied  to  the  surface  dQ.  at  position  rs  and  S  is 
expressed  as: 

S  -  S8(F-rJ  (14) 


where  rs  is  the  position  vector  of  the  point  source.  By  definition,  we  have 


(15) 


{  S8(r  -  rs]dT  =  (•?  r'  *r* 
rei  [0  otherwise 

where  el  is  the  element  of  the  finite  element  mesh.  From  Eqn(12)  we  have 

rel  =-  I  LjS5(?-^)ir  =— Lj(rs)S 
Y  pel  Y 

If  the  position  rs  of  the  point  source  coincides  with  the  position  of  node  j  ,  when  Lj(rs )  =  1 , 
i.e.  we  have 


2.3  Gale  rkin  formulation  of  finite  element  solution  of  &m 

The  forward  finite  element  solution  of  the  emission  diffusion  equation  is  formulated  similarly 
to  the  excitation  diffusion  equation  with  the  exception  that  (i)  the  complex  excitation  fluence 
must  be  first  computed,  and  that  (ii)  the  formulation  of  fluorescence  lifetime  poses  some 
additional  considerations.  As  in  the  formulation  of  the  solution  for  Ox,  parameters  of  , 

T,  and  Dm  are  expected  to  vary  linearly  within  the  triangular  elements.  In  order  to  facilitate 
the  formulation  of  the  Galerkin  finite  element  method,  the  RHS  term  of  Eqn  (2)  is  expressed 
in  terms  of  a  binomial  expansion  assuming  m  <1  and  higher  order  terms  are  neglected. 
Hence,  the  RHS  source  term  of  Eqn  (2)  becomes  Sm  [l  +  /co(r)t]<& v(r,co). 

This  assumption  is  appropriate,  since  the  phase-delay,  0m,  of  the  generated  fluorescent 
intensity  modulated  wave  reaches  a  maximum  of  90  degrees  with  respect  the  local  excitation 
intensity  modulated  wave  and  becomes  insensitive  to  fluorescence  lifetime  when  coot.  The 
formulation  of  the  emission  finite  element  solution  for  is  identical  to  that  presented  in 

section  2.2  for  Ox  with  the  exception  that  the  parameters  Om,  Dnl,  ,  and  x  are 

assumed  to  vary  linearly  within  each  element.  The  finite  element  formulation  of  emission 
diffusion  is  given  below: 


(16) 
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It  is  important  to  emphasize  that  the  absorption  owing  to  fluorophore  is  equivalent  in  our  case 
to  \in  and  to  \in  .  We  differentiate  between  the  two  parameters  when  reconstructing 

r  Ux_^m  r  w  .r  f 

absorption  from  synthetic  excitation,  ,  and  absorption  from  emission  data  |l^_^  . 


2.4  Solution  of  the  system  of  finite  element  equations 

For  both  excitation  and  emission  finite  element  formulations,  the  local  stiffness  matrices, 
K/f,  are  formed  and  then  combined  into  a  complex,  global  stiffness  matrix,  K,  which  depends 
upon  the  absorption  and  scattering  properties  at  each  wavelength.  For  solution  of  excitation 
fluence,  <PX ,  the  vector  b  reflects  the  complex  source  of  excitation  light  Irom  boundary  point 

sources  at  positions  rs,  while  for  solution  of  the  emission  fluence,  0/;, ,  vector  b  represents  the 
source  generated  from  within  volume 

To  solve  for  O*  m ,  the  complex  set  of  sparse,  linear  equations  generated  by  the 
finite  element  method, 


K*,.«  =  b  (17) 

are  efficiently  solved  using  a  subroutines  ZGBTRF  and  ZGBTRS  (LAPACK  subroutines). 
ZGBTRF  computes  a  LU  factorization  of  a  complex  M2  by  N  band  matrix  using  partial 
pivoting  with  row  interchanges.  (M2  =  2*KL+KU  +  1 ,  KL  and  KU  are  the  numbers  of 
subdiagonals  and  superdiagonals  within  the  band  of  K,  respectively).  ZGBTRS  solves  a 
system  of  linear  equations  with  a  general  band  matrix  K  using  the  LU  factorization  computed 
by  ZGBTRF.  The  main  advantage  of  this  method  is  that  once  the  decomposition  is  done,  the 
solution  vector  can  be  obtained  for  any  number  of  right  hand  side  vectors.  The  accuracy  of 
the  solution  for  Ov<m  depends  crucially  on  the  mesh  size,  A,  with  the  accuracy  improving  with 

h\ 


з. 0  Inverse  solution  for  fluorescence  FDPM 

The  variables  in  the  inverse  problems  are  the  absorption  coefficients,  (for  excitation), 

и .  (for  emission),  and  the  lifetime  T  at  each  nodal  point  of  a  finite  element  model.  If  there 

'  ax—>m  v 

are  N  nodal  points  then  the  number  of  measurements  must  be  greater  than  N.  The  number  of 
measurements  made  at  N  B  boundary  nodes  for  Ns  sources,  requires  that 


(NB-l)Ns>N  (IB) 

Due  to  the  reciprocality  of  sources  and  detectors,  i.e.,  when  a  detector  position  is  swapped 
with  a  source  position  equivalent  measurements  are  registered,  then  we  require: 

(jVg  -  >  N  (19) 

2 

If  /  =  l,...,iVJ  denotes  the  different  distribution  of  source  s  and  j  =  denotes  the 

boundary  nodes  at  which  measurements  are  made,  the  error  function  for  optimization  of 
absorption  and  lifetime  imaging  may  be  taken  as 


(20) 


The  subscript  c  denotes  the  values  calculated  by  the  forward  simulator  problem  and  the 
subscript  me  denotes  the  experimental  (or  in  this  case  synthetically  generated)  fluence  values. 
The  superscript  *  denotes  the  complex  conjugate  of  the  complex  function  The 

subscript  /  and;  denote  fluence  values  that  arise  between  source  /  and  detector;'. 

It  is  seen  from  the  equation  (20)  that  the  calculation  of  error  function  E  involves  N s 

solutions  of  the  direct  problem.  The  global  stiffness  matrix  K  v  m  for  excitation  and  emission 

remains  constant  with  only  the  vector  b  changing  with  each  excitation  source.  Since  the 
global  stiffness  matrix  Kx  depends  on  jla^  and  b  is  independent  of  j ia^,  LU 

decomposition  of  Kx  for  solution  of  <E>_V  need  only  be  performed  once  in  each  iteration  of 
the  optimization  problem. 

An  efficient  way  of  calculating  Ex  m  for  excitation  and  emission  equations  consists 
of  the  following  steps: 

Step  1  Given  \ia  ,  calculate  the  local  stiffness  matrices  KW(jI^  ) 


Step  2 
Step  3 

For  each  source  / 
Step  4 

Step  5 


Step  6  £x>m  =  -  X 


K=£Kel  and  set  F  =  0 
el 

Decompose  K  into  the  LU  decomposition  LU  =  K 


calculate  b(s) 

Solve,  LU$m  =b  gives 
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As  most  of  the  efficient  optimization  algorithms,  the  truncated  Newton  method  requires  the 
function  value  Ex  m  and  the  gradient  of  the  error  function,  V£r  m ,  in  order  to  efficiently 

update  parameter  estimates.  The  gradients  of  the  error  function  in  Eqn  (20)  represent  the 
direction  in  which  the  error  function  is  increasing/decreasing  most  rapidly.  Since  absorption 
coefficients  ,  ]ln  and  Tare  the  unknown  parameters,  the  gradients  are  obtained  by 

r  axf  r  ux->m 

taking  partial  derivatives  of  Eqn  (20)  with  respect  to  \ia^  ,  (Xar_>m  and  t. 


( 


fk) 

- 

c 

k) 

1 

me 

me  j 

’  / 

1  / 

(21) 


denotes  the  inner  product  of  two  complex  vectors  and  Re(.)  denotes  the  real  part  of  a 

complex,  number.  The  gradient  of  the  real  valued  error  function  with  respect  to  real  valued 
rr  II  ,  and  x  remains  real  valued.  At  minimum  these  gradients  should  be  zero,  so  the 

real  part  of  the  gradients  are  zero.  Therefore,  we  consider  only  the  real  part  of  the  gradients. 
In  the  following,  the  optimization  strategy  using  the  truncated  Newton  method  with  the 
gradients  computed  by  reverse  automatic  differentiation  is  described. 

3.1  Truncated  Newton's  method  as  an  optimization  strategy  for  FDPM  imaging 

The  truncated  Newton  method  with  trust  region  is  used  for  non-linear  optimization  problem. 
The  reader  interested  in  background  literature  is  referred  to  the  references  [27,  28]. 

Newton’s  method  is  based  upon  approximating  the  function  Exm  (jv  +d) 

( m  denotes  p.a  {for  excitation}  or  (for  emission}  in  this  section)  at  the  Mi  iteration 

by  the  quadratic  model: 

+  (H-o*  )  +  Q(^) 


2(d)  =  g[d  +  ^drG*d  (24) 

In  these  equations  g*=g (jXa;.  )=V£^m(jIaj J  and  G*  =V2£Y  )  denote  the  gradient 

vector  and  Hessian  matrix,  respectively.  The  Newton  direction  is  obtained  from  an  exact 
minimum  of  above  equation;  i.e.  the  search  direction  d  at  iteration  k  is  defined  by  the 

equation 

G*d  =  -g*  (25) 

A  sequence  of  approximate  solutions  of  the  above  equation  is  generated  by  the  conjugate 
gradient  method  until  a  vector  d,  (see  Step  3)  is  obtained  for  which  the  following  condition 
on  the  relative  residual  is  satisfied  as  suggested  by  Dembo  and  Steihaug  [29]; 


where  r =  G^d f-  +  gk  .  The  conjugate  gradient  iteration  is  then  truncated  and  d,  is  used  as 

the  search  direction  for  the  minimization.  The  relative  residual  is  used  as  a  measure  of 
accuracy  and  the  required  level  of  accuracy  improves  as  the  minimization  proceeds  and  as 
(jl^  j  approaches  a  minimum.  It  is  not  necessary  to  compute  the  Hessian  G  as  only  the 

product  Gd  is  required.  The  product  is  calculated  using  the  finite  difference  formula  given 
below: 

G(x)d  =  —  [g(x  +  ad)  -  g(x)]  (26) 

a 


where  a  = 


j  machine  precision 


.  This  avoids  the  calculation  and  storage  of  the  Hessian  but 


requires  an  additional  gradient  evaluation  for  each  minor  iteration.  The  pseudo  codes  for 
truncated  Newton  optimization  is  as  follows: 


Step  1  Initialization 

An  initial  guess  (j^  )  is  made  of  the  solution:  The  radius  R{  =  0.01  of  the  trust 
region  is  defined  around  £0“  £(p^0)  g0  =  are  comPuted- 

Step  2  Stopping  Criteria 

If  ||g;.||  <  e  exit,  e  is  specified  by  the  user 

Else  continue  with  next  Step  3 
Step  3  Computing  Newton  Direction 

Apply  conjugate  gradient  method  to  find  an  approximate  solution  d*  of  the  Newton 
equation  Gkdk  =  -gk  such  that  |djt||  <  Rk  and  d[%k  <  -ei[|d*||||g*[|>  typically 
ei  =  0.001  [30].  The  conjugate  gradient  algorithm  ensures  that  at  the  y-th  iteration 
d j+i  =d j  +  a jp j  remains  within  the  trust  region  (a }  =  steplength  and  p }  is  the 

conjugate  direction) 

Step  4  Line  Search 

After  computing  the  truncated  Newton  direction  d^  ,  a  line  search  is  used  to  find  an 
appropriate  Xk  that  reduces  the  function  Ex  m  along  the  line  (1^  +kkdk  .  The 

Armijo  line  search  [31]  is  used.  It  calculates  Xk  such  that  Wolfe’s  [32]  conditions  2 

and  3  are  satisfied. 

Condition  2: 

Efca  +  ad)  -  £’(pa)  ^  e2cxdrg  where  Zo  ls  a  constant,  such  that 
0.5>82>0  (typically  82=0.1). 

Condition  3 

E(jia  +2ad)  >  +2e3adrg)  where  e3  is  a  constant  such  that 
0.5  >  e3  >  0 

Step  5  New  Approximation 

Compute  Ffljt+l  =  \xak  +  Xkdk  ,  £x?%+1 ,  g*+i 

Step  6  Update  the  trust  region 

Once  Xk  has  been  found  the  radius  of  the  trust  region  Rk+[  is  altered  as  follows: 


fl,=0.01 

&k+l  =  2 Rk  if  ^  *  10 
^Jt+l  =  -J**  if  *•*  <1-0 

The  truncated  Newton  search  direction  is  constructed  so  that  it  always  satisfies  Wolfe’s 
condition  1  (Step  3,  dTkgk  <  -ei||d*||||g*|| ).  Wolfe’s  conditions  2  and  3  are  satisfied  within 

the  line  search,  as  discussed  is  Step  4.  Since  all  the  Wolfe’s  conditions  have  been  satisfied, 
the  method  is  globally  convergent,  and  will  terminate  in  a  finite  number  of  iterations. 

A  similar  formula  applies  for  lifetime  x  and  absorption  coefficient  |Aa^w  . 


3.2  Reverse  automatic  differentiation  for  computation  of  V£xm 


The  reverse  automatic  differentiation  (RAD)  method  is  used  to  calculate  the  gradient 
g  =  g(jlfl  or  x  or  .  The  reader  interested  in  background  literature  is  referred  to 

References  [28,  33,  and  34].  A  brief  description  of  the  reverse  differentiation  and  analytic 
implementation  of  this  method  is  discussed  in  the  following  section. 

Reverse  automatic  differentiation  is  discussed  in  terms  of  the  Wengert  [35]  list.  This 
list  decomposes  complicated  function  of  many  variables  into  a  sequence  of  simpler  functions 
of  one  or  two  variables.  Functions  of  two  variables  are  called  the  "binary  function".  Binary 
functions  are  either  addition,  or  subtraction  or  multiplication.  Functions  of  one  variable  are 
either  reciprocation,  raise  of  power,  exponential,  logarithm,  trigonometric,  or  hyperbolic. 
These  functions  are  defined  as  unary  functions. 

If  f(x)  is  a  function  of  the  n  variables  jq, then  a  set  of  new  variables 

are  introduced,  where  P  is  the  number  of  arithmetic .  operations  involved  in 
calculating  the  function  f(x).  A  Wengert  list  for  the  calculation  for  the  calculation  of  f(x) 
consists  of  a  list  of  these  unary  and  binary  function  processed  in  a  given  order; 


Given 

For 

then  if  Ft  is  binary 
and  if  F,  is  unary 


Jt,-,  /  = 

i  =  n  +  l,...,F 
*«  =  Ft  (Xj ,  xk  )  j,  k  <  i 
Xi=F{xjl.  j  <  i 
fix)  =  xn+P 


Reverse  differentiation  can  be  explained  as  follows.  Given  the  Wengert  list  for  the  calculation 
off  a  set  of  adjoint  variables  xt  =  is  introduced  and  the  order  of  the  operations 

reversed.  The  method  is  based  on  the  function  of  a  function  rule  of  calculus  that  implies 


;>f 

dx{  j  dxj  dxj 


The  automatic  method  involves  a  pass  in  the  reverse  direction  through  the  list  after  the 
function  has  been  evaluated  in  a  forward  sweep.  Hence  the  method  is: 


Given  x,- ,  i  =  I 

set  .?/  =0,  i  =  1,...,«  +  P-1 

for  /  =  n  +  P,...,«  +  l 

then,  if  Z7*  is  binary,  Xj  =  Xj 

and  xk  =  xjf, 

else  if  Z7  is  unary,  xj  =  x; 

derivatives  g,-  =  x/ 


and  xn+P 

dFt 


+  Xi 


dX: 
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dxk 
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■ 1  a,, 

i  =  l,...,/i 


=  1 

i  >  J 
i  >  J 
i  >  j 


Griewank  [14]  showed  that  it  is  always  possible  to  calculate  any  gradient  vector  of  a  function 
in  less  than  3  times  the  computational  cost  of  the  function  by  the  reverse  automatic 
differentiation  method  irrespective  of  the  dimension  of  the  problem.  The  main  disadvantage 
of  the  automatic  code  is  that  it  requires  large  storage  and  significant  data  accessing  overheads 
[28,  33,  and  34].  In  order  to  overcome  these  difficulties,  the  reverse  differentiation  method  is 
extended  and  has  been  performed  analytically.  This  does  not  require  large  storage  and  at  the 
same  time  gradients  are  calculated  in  less  than  three  times  the  operation  count  of  the  function. 


The  gradients  of  Ex  m  with  respect  to  parameter  \ia ^  are  calculated  using  the  reverse 
differentiation  method  as  follows. 


Step  1 
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Step  2  Calculate  element  matrix  KeI  (jia ) 

Step  3  Assemble  global  matrix  K 

Step  4  Decompose  K  =  LU 

For  each  light  source  s; 

Step  5  Calculate  b(s) 


Step  6 
Step  7 
Step  8 

Step  9 


Solve  LUO*,,,  =  b(s)  to  give  (®.v,m)c 
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for  all  j<£  d£l 


Step  10  Solve  L\Jvx  m  =  <Px  t] 


Step  11  Update  K  =  K  -  vT x,m®x ,m,  b  =  b  + 


Step  12 


The  above  gradients  of  error  function  £  Y  mare  calculated  using  the  global  stiffness  matrices, 

K,  and  the  matrix  b  assembled  from  the  excitation  and  emission  diffusion  equations  as 
described  in  Section  2.2  and  2.3.  (The  reverse  differentiation  method  described  above  uses  the 

fact  that  if  0Am  is  the  solution  of  the  equation  KOY<m=bthen  K  =-  vr.Y,mOv  m  and 

b  =  vvm  where  vxm  is  the  solution  of  the  equation  KvAV„=  Oy/„  ,  see  Appendix  A).  With 
VEy  m  and  Ex  m  computed,  the  Newton  directions  can  be  computed  as  described  in  Section 
3.1.  Similar  formulas  hold  for  x  and  jL 

4  Conclusions 

We  have  formulated  the  solution  to  the  absorption  and  fluorescence  FDPM  optical  imaging 
problem  as  a  non-linear  optimization  scheme.  We  develop  the  truncated  Newton  method  for 
minimizing  the  sum  of  errors  between  measured  and  synthetic  data  generated  using  a 
simplistic  finite  element  model  with  a  minimal  number  of  sources  and  detectors.  Since  the 
truncated  Newton  method  requires  only  the  product  of  the  Hessian  with  a  direction  at  any 
point,  it  is  computationally  efficient  for  large  optimization  problems.  We  couple  the  truncated 
Newton  method  with  a  finite  element  solver  which  assumes  a  polynomial  series  in  absorption 
and  lifetime  in  order  to  accurately  compute  Oxm  and  employ  principles  of  reverse 
differentiation  in  order  to  efficiently  compute  the  gradient  of  the  error  function,  VEX  m  . 

Using  synthetic  data  with  simulated  measurement  noise,  we  demonstrate  in  Part  II  the 
optimization  strategy  on  a  series  of  studies  involving  a  minimum  number  of  sources  and 
detectors. 


Appendix  A 

Calculation  of  adjoint  variables 

In  this  appendix  we  derive  the  expressions  for  the  adjoint  variables  = 
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(since  K  is  a  function  of  \ia  ,  we  differentiate  K  with  respect  of  \La ^  ).  Now  we  use  the 
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Since  the  source  term  b  is  constant  in  excitation  equation. 
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On  substituting  Eqn  (A2)  in  Eqn  (Al)  we  have 


(A2) 
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where 


v  =  oat1 
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We  can  solve  the  Eqn(A4)  and  obtain  v  .  Finally  K  is  obtained  using  Eqn(A3). 
(ii)  Adjoint  variable  b 

3£,.m  _3£,.r.m  =  £  5g>V.m  _  ±  3<t>.v.m  9X 


b  = 


-=  0 r 


at  a*  at 


at 


at  at 


(A5) 


(since  t  is  a  function  of  t ,  we  differentiate  b  with  respect  of  T ).  Now  we  use  the  global 
—  9x 

equation  K<t>x  m  =b  to  find  —  as  follows: 
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Substitute  Eqn(A6)  into  Eqn(A5)  we  have 
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Eqn(B7)  is  the  same  as  Eqn  (B4).  Hence  b  =  v 

Nomenclature 


M  amplitude  modulation 

r  position 

c  speed  of  light 

S(rs  co )  source  at  location  rs  and  modulated  frequency  co . 

D  optical  diffusion  coefficient  -  1/3  [jla  +|ij 

n  normal 

rs  source  position 

rd  parameter  describing  y 

nrel  relative  refractive  functions 

Wj  independent  weighting  functions 

Lj  co-ordinates  for  triangular  elements 

K?  element  stiffness  matrices 

rel  local  vector  in  FEM  stiffness  equations 

b  local  vector  in  FEM  stiffness  equation  containing  source  terms 

Z  variance 

G(0,1)  random  number  with  Gaussian  distribution  of  zero  mean  and  unit  variance 

N  number  of  nodal  points 

Ns  number  of  sources 

N B  number  of  boundary  nodes 

p  distance  between  source  and  detector  pair 

0  phase-delay 

co  angular  modulation  frequency 

|ia  absorption  coefficient 

jLig  isotropic  scattering  coefficient 

<I>  complex  fluence 

<j>  quantum  efficiency 

t  fluorescence  lifetime 

H  volume 

y  parameter  to  account  for  refractive  index  mismatch 

T  surface 

8  stopping  criteria  for  truncated  Newton’s  method 

E  error  function  for  all  sources 

F  error  function  for  individual  sources 

gradient  vector  at  iteration  k  =  VEk 


d  search  direction 

Gk  Hessian  matrix 
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Abstract:  Using  two  dimensional  synthetic  frequency-domain  measurements,  the 
inverse  imaging  problem  is  solved  for  absorption  and  fluorescence  lifetime  mapping 
with  the  truncated  Newton’s  optimization  scheme  developed  in  Part  I  of  this 
contribution.  Herein,  we  present  reconstructed  maps  of  absorption  owing  to  a 
fluorophore  from  excitation  and  emission  measurements  which  detail  the  presence  of 
tissue  heterogeneities  characterized  by  tenfold  increase  in  fluorescent  contrast  agent. 
Our  results  confirm  that  fluorescence  provides  superior  mapping  of  heterogeneities 
over  excitation  measurements.  Using  emission  measurements  we  then  map 
fluorescent  lifetime  under  conditions  of  tenfold  uptake  of  contrast  agent  in  tissue 
heterogeneities.  The  ability  to  map  fluorescent  quenching  and  lengthening  of 
contrast  agents  facilitates  the  solution  of  the  inverse  problem  and  further  improves 
the  ability  to  reconstruct  tissue  heterogeneities. 
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1.  Introduction 

In  Part  I  of  this  contribution,  we  have  formulated  the  absorption  and  fluorescence  inverse 
imaging  problem  as  an  optimization  problem  in  which  an  interior  map  of  tissue  optical 
properties  of  absorption  and  fluorescence  lifetime  can  be  reconstructed  from  excitation  and 
emission  frequency-domain  photon  migration  (FDPM)  measurements.  Since  significant 
challenges  for  the  non-linear  and  large  scale  reconstructions  are  further  exaggerated  by  the 
physiological  low  contrast  of  absorption  and  scattering  property  differences  between  normal 
and  diseased  tissues  [2],  we  have  focused  upon  the  development  of  absorption  and 
fluorescence  lifetime  imaging  following  the  administration  of  a  fluorescent  contrast  agent. 
Previously,  we  have  conducted  FDPM  experimental  measurements  to  show  that  the  FDPM 
contrast  owing  to  fluorescence  exceeds  that  possible  by  monitoring  absorption  with  excitation 
FDPM  measurements  [3].  One  objective  of  this  contribution  is  to  compare  the  reconstruction 
of  absorption  cross  section  owing  to  a  fluorescent  agent,  which  has  a  tenfold  preferential 
uptake  in  the  tissue  region  of  interest  using  excitation  and  emission  FDPM  measurements. 
Secondly,  we  focus  upon  the  reconstruction  of  fluorescent  lifetime  for  cases  in  which  a 
tenfold  preferential  uptake  exists  within  the  tissue  volume  of  interest  is  accompanied  by  a 
shortening  and  lengthening  of  a  first  order  fluorescence  decay  time  (or  fluorescence  lifetime). 
In  the  following,  we  describe  the  generation  of  the  synthetic  data  set  for  input  into  the 
optimization  scheme,  and  then  the  results  for  absorption  and  fluorescence  lifetime  imaging. 
Finally,  we  comment  on  the  availability  of  fluorescent  contrast  agents  for  inducing  optical 
contrast  and  the  future  work  to  adapt  our  inversion  strategies  to  actual  FDPM  measurements. 

2.0  Generation  of  synthetic  data  set  using  forward  simulator 

In  order  to  test  our  inversion  a  truncated  Newton’s  optimization  scheme  developed  in  ref.  [1], 
we  generated  synthetic  data  of  O x  m  in  a  two  dimensional,  square  domain  of  4  x  4  cm2.  For 

the  forward  simulation,  the  mesh  contained  2209  internal  and  192  boundary  nodes  with  4608 
triangular  elements.  It  should  be  noted  that  a  uniform  grid  is  used  in  this  exercise  whereas  for 
future  work  we  adapt  a  finer  grid  near  the  source  points  where  the  gradients  are  very  high. 
Four  sources  of  intensity  modulated  excitation  light  were  simulated  midpoint  on  each  side  of 
the  domain  with  a  total  of  59  point  detectors  simulated  equidistant  along  the  domain 
periphery  excluding  positions  occupied  by  one  of  the  four  sources  and  at  each  corner  of  the 
domain.  Predictions  of  <&x  m  were  made  at  each  detector  for  each  of  the  four  sources,  each 

providing  intensity  modulation  at  50,  100  and  150  MHz  with  unit  depth  of  modulation  and 
zero  phase.  For  this  simulated  measurement  configuration  there  were  4  x  59  =  236  simulated 
observations  obtained  for  reconstruction  of  interior  optical  properties  at  the  excitation  and 
emission  wavelengths.  At  the  excitation  wavelength,  “heterogeneities”  were  to  be  detected 
based  upon  their  absorption,  and  at  the  emission  wavelength,  based  upon  their  fluorescence 
lifetime.  Table  I  lists  the  background  and  “heterogeneity”  optical  properties,  as  well  as  the 
location  and  size  of  the  simulated  “heterogeneities”  used  for  the  generation  of  synthetic  data 


for  reconstruction.  The  simulator  was  coded  in  Fortran  77  and  took  6  seconds  on  a  SUN 
Ultrasparc  10  Workstation  (200MHz). 


Table  1.  Parameters  used  in  truncated  Newton  method  optimization  for  absorption  and  lifetime  imaging 


background 

targets 

reconstruction 

variables 

cm”1 

or 

cm-1 

cm”1 

cm-1 

X 

ns 

0 

or 

cm”1 

X 

ns 

case  1, 

Figure  1 

0.0 

0.02 

10.0 

N/A 

N/A 

N/A 

0.2 

N/A 

case  II,  Li 

Figure  2 

0.0 

0.02 

10.0 

0.02 

1.0 

0.034 

0.2 

10.0 

case  III,  X 

Figures  3  &  4 

0.0 

0.02  | 

10.0 

0.02 

1.0 

0.034 

0.2 

10.0 

case  IV,  X 

Figures  5  &  6 

0.0 

0.02 

10.0 

1.0 

Footnotes:  case  I-  excitation  equation  is  used,  emission  equation  is  used  for  all  other  cases.  In  this  study,  we  assume 
all  absorption  at  the  excitation  wavelengths  is  owing  to  fluorophore. 

To  mimic  measurement  error,  zero-mean,  Gaussian  noise  corresponding  to  1.0% 
standard  error  in  fluence  was  added  to  the  synthetic  data  sets  by  using  the  formula: 

d>r  =  $'(l.0  +  Z*G(0,l))  (1) 

where  G(0,  l)  is  a  Gaussian  distribution  with  zero  mean  and  unit  variance;  Z  =  0.01 ;  and 
O'  are  the  simulated  data  without  noise. 


To  simulate  0.1  degree  noise  in  phase,  we  compute  the  error  in  the  final  complex 
fluence,  <&'x : 


RI 

_ k) 

1  -  tan9tan(0.1*G)  re  O1,.) 


tan(0  +  0.1*G(0,l))  = 
tan  0  +  tan(0.1*G) 


img 

re 

img 


+  tan(0.1*G)  ,  , 

rej^x)  _  _  j 

|  -  !&&jd  ,lan(0.^  '  /»&) 

n?($x) 

fmg(<&Y)  +  re(Ox)*tan(0.1*G)  _  /mg(ojY ) 
re(O^)  -  img (O  *  )*  tan  (0. 1  *  G ) 
imgfpx)**  img(®x)  +  re(<f>JC)*tan(0.1*G) 
re(oJY)  =  re(<&x)  -  /mg(Ov)*  tan(0.1*G) 


Now  our  new  fluence  is:  &'x  =  (re(o^. )  img  |ct>  Y  ])  (2) 

The  synthetic  fluence  for  emission  is  similarly  generated  and  employed  as  input  to  the 
iterative  inversion  algorithm. 

A  user  specified  double  precision  parameter  £  is  used  as  the  convergence  parameter  within 
the  inversion  scheme  (see  section  3.1  of  reference  [1]).  The  iterative  method  is  stopped  if  the 
length  of  the  gradient  vector  is  less  than  e .  The  final  images  reported  herein  are  the  results  of 

iterations  until  the  length  of  the  gradient  ||g||  is  less  than  £  =  1CT9  ,  chosen  by  trial  and  error. 

3.0  Results  and  discussion 

A  coarse  mesh  with  225  internal  nodes  and  64  boundary  nodes  (total  289  nodes)  is  used  for 
the  inversion  problems,  whereas  total  2401  nodes  with  uniform  mesh  is  used  for  forward 
problem.  Using  the  synthetic  data  sets  described  in  Section  2,  we  demonstrate  the  truncated 
Newton  method  for  reconstruction  of  absorption  from  excitation  and  fluorescence  FDPM 
measurements  as  well  as  for  reconstruction  of  lifetime  from  fluorescence  FDPM 
measurements.  For  this  mesh  one  function  (Ex^m ) evaluation  took  9.1  seconds  and  one 
function  and  gradient  calculation  by  reverse  differentiation  took  10.5  seconds. 

3.1  Absorption  imaging  from  synthetic  excitation  measurements 

As  described  in  Table  I,  Case  I  synthetic  data  set  consisted  of  detecting  three 
“heterogeneities”  based  upon  a  ten-fold  increase  in  absorption  as  might  occur  upon  the 
administration  of  a  contrast  agent.  Figure  1(a)  illustrates  the  actual  distribution  of  absorption 
within  the  phantom  with  a  background  absorption  \\.a^  value  of  0.02  cm’1.  Upon  using  the 

background  absorption  value  as  an  initial  starting  guess,  our  attempt  to  reconstruct  the  spatial 
distribution  of  absorption  with  simulated  measurement  with  noise  was  aborted  after  25 
iterations  because  the  results  remain  unchanged.  The  result,  illustrated  in  Figure  1(b),  shows 
that  the  heterogeneities  may  be  located  and  differentiated  from  the  relatively  uniform 
background.  However,  the  differentiation  of  the  heterogeneities  from  one  another  and  the 
quantitation  of  their  ten-fold  increase  in  absorption  are  disappointing.  The  average  values  of 
the  reconstructed  absorption  estimate  of  each  of  the  three  “heterogeneities”  are  plotted  as  a 
function  of  iteration  in  Figure  1(c).  While  oscillations  in  the  parameter  estimates  occur  in  the 
initial  iterations,  the  average  values  appear  to  smoothly  approach  an  absorption  value  that 
underestimates  their  actual  values.  These  numerical  approximation  errors  may  be  due  to 
inadequate  discretization  of  the  model  geometry  by  the  finite  element  method.  Since  the 
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Figure  2  (a)  Reconstructed  absorption,  |ifl  ,  from  fluorescence  measurements 

and  (b)  Average  value  of  | |Xfl  .as  a  function  of  iteration. 


the  synthetic  fluorescence  FDPM  measurements  as  described  by  case  II  in  Table  I.  Figure 
2(b)  shows  the  corresponding  average  values  within  each  of  the  heterogeneities  as  a  function 
of  iteration.  As  expected  from  the  superior  contrast  offered  by  fluorescence  over  absorption 
and  the  linearity  of  emission  equation  (see  Eqn  (2)  of  ref  [1]),  the  convergence  upon  an 
absorption  map  from  fluorescence  measurements  occurs  faster  and  more  accurately  than  from 
excitation  measurements.  Reconstruction  of  absorption  from  excitation  requires  recalculation 
of  the  global  matrix  at  each  iteration  while  the  global  matrix  in  absorption  or  lifetime 
reconstructions  from  emission  remains  unchanged  throughout  the  calculation. 

3.3  Lifetime  imaging  from  synthetic  fluorescence  FDPM  measurements 

While  the  contrast  owing  to  absorption  may  provide  localization  of  tissue  disease,  the 
discrimination  of  diseased  tissues  based  upon  changes  in  lifetime  of  exogenous  fluorescent 
contrast  agents  has  been  demonstrated  in  endoscopic  applications  [4]  and  proposed  for  optical 
tomography  [5-7].  In  the  following  two  studies,  we  demonstrate  the  truncated  Newton  method 
for  demonstrative  cases  of  fluorophore  quenching  (shortening  of  fluorophore  lifetime)  and 


enhanced  activated  fluorophore  stability  (i.e..  lengthening  of  fluorophore  lifetime)  within 
tissue  “heterogeneities.” 

3.3.1  Imaging  based  upon  lengthening  of  fluorophore  lifetime 

Figure  3(a)  depicts  the  actual  lifetime  map  of  the  Case  III  study  listed  in  Table  I  in  which 
three  heterogeneities  with  ten-fold  uptake  of  fluorescent  dye  exhibit  a  lengthening  of  lifetime 
(xb=10  ns)  within  a  background  which  possessed  a  shorter  lifetime  (th=l  ns).  With  an  initial 
starting  guess  equal  to  the  background  lifetime,  the  reconstructed  maps  of  lifetime  are  shown 
for  50MHz  (Figure  3(b));  100  MHz  (Figure  3(c));  and  150  MHz  (Figure  3(d)).  The  quality  of 
the  reconstruction  map  seems  to  improve  with  increasing  modulation  frequency  and  Figure 
4(a)  through  (c)  suggests  more  rapid  convergence  at  increased  modulation  frequencies.  It  is 
noteworthy  that  using  the  Levenberg-Marquardt  approaches.  Paithankar,  et  al  [5]  and  Jiang 
[8]  report  only  reconstruction  involving  fluorophore  quenching  and  were  unable  to 
successfully  reconstruct  heterogeneities  with  longer  lifetimes  within  a  background  of  short 
lifetime.  Two  additional  observations  of  the  results  can  be  made:(i)  while  the  recovered 
values  of  the  absorption  from  excitation  measurements  underpredict  their  local  absorption 
coefficient  (0.05  as  opposed  to  0.20  cm'1),  the  recovered  values  of  lifetime  appear  to  be  closer 
to  the  actual  values  (7  instead  of  10  ns);  (ii)  while  the  absorption  imaging  from  excitation 
FDPM  measurements  (Figure  1)  did  not  achieve  convergence  after  25  iterations,  convergence 
of  fluorescence  lifetime  imaging  was  achieved  in  20-30  iterations.  As  described  below  similar 
results  were  obtained  for  fluorescence  lifetime  imaging  involving  fluorophore  quenching. 


(c)  (d) 


Figure  3  (a)  "True'  distribution  of  fluorophore  lifetime,  possessing  longer  lifetime  within 
three  heterogeneities  having  ten-fold  uptake  of  fluorescent  dye;  (b)  Reconstructed  lifetime 
at  50  MHz;  (c)  100  MHz;  (d)  150  MHz. 


3.3.2  Imaging  based  upon  fluorophore  quenching 


Figure  5(a)  depicts  the  actual  lifetime  map  of  the  Case  IV  study  listed  in  Table  I  in  which 
three  heterogeneities  with  ten-fold  uptake  of  fluorescent  dye  exhibit  a  shortening  of  lifetime 
(xb=l  ns)  within  a  background  which  possessed  a  longer  lifetime  (xh=10  ns).  We  assume  all 
other  parameters  are  known.  With  an  initial  starting  guess  equal  to  the  background  lifetime, 
the  reconstructed  maps  of  lifetime  are  shown  for  50MHz  (Figure  5(b));  100  MHz  (Figure 
5(c));  and  150  MHz  (Figure  5(d)).  The  quality  of  the  reconstruction  map  again  seems  to 
improve  with  increasing  modulation  frequency  and  Figures  6(a)  through  (c)  again  suggests 
rapid  convergence  at  increased  modulation  frequencies  and  greater  efficiency  than  can  be 
shown  with  absorption  imaging. 


(c) 
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Figure  5  (a)  'True'  distribution  of  fluorophore  lifetime,  with  fluorophore  quenching  within 
three  heterogeneities  having  ten-fold  uptake  of  fluorescent  dye;  (b)  Reconstructed  lifetime 
at  50  MHz;  (c)  100MHz;  (d)  150  MHz. 


Reconstruction  of  lifetime,  background  1  ns,  targets  10  ns,  50Mhz 


Reconstruction  of  lifetime,  background  1  ns,  targets  10  ns,  lOOmhz 


(b) 


Reconstruction  of  lifetime,  background  1  ns,  target  ns,  150  mhz 


number  of  iteations 

Figure  6  The  distribution  of  fluorophore  lifetime,  quenching  fluorescent  lifetime  within  three 
hetergeneities  having  ten-fold  uptake  of  fluorescent  dye.  Average  value  of  heterogeneity 
lifetime  as  function  of  iteration;  (b)  at  50  MHz;  (c)  100  MHz;  (d)  150  MHz 
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Conclusions  and  Future  Work 


In  this  contribution,  we  demonstrate  the  truncated  Newton's  optimization  scheme  for 
absorption  and  fluorescence  optical  tomography.  We  have  presented  a  method  to  calculate  the 
gradient  of  error  function  based  on  reverse  differentiation  method  in  a  finite  element  method 
based  scheme.  This  method  is  programmed  by  hand  so  that  the  overhead  problems  usually 
associated  with  automatic  differentiation  do  not  occur.  It  is  shown  that  a  function  and  the 
gradient  calculations  are  cheaper  than  three  function  evaluations  using  reverse  differentiation. 
Our  work  confirms  that  reconstruction  of  absorption  owing  to  a  contrast  agent  is  enhanced  if 
fluorescence  measurements  are  made  over  excitation  measurements.  In  addition,  the  ability  to 
detect  heterogeneities  with  a  tenfold  uptake  of  dye  that  experiences  lengthening  and 
shortening  of  fluorescence  lifetime  is  demonstrated.  In  practice,  FDPM  optical  tomography 
must  be  accomplished  with  a  substantive  number  of  source-detector  measurements  and  for 
three-dimensional  geometries  [9,10].  While  work  continues  to  develop  three-dimensional 
multi-pixel  FDPM  measurements  [11]  using  fluorescent  contrast  agents  [12],  the  development 
of  effective  inversion  strategies  to  handle  large  data  sets  at  minimal  computational  cost  and 
storage  burden  is  paramount.  We  believe  the  continued  development  of  gradient-based 
optimization  schemes,  such  as  presented  herein,  are  necessary  for  stable  recovery  of  interior 
optical  property  maps  for  medical  imaging. 
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ABSTRACT 

Near-infrared  biomedical  optical  imaging  consists  of  imaging  interior  volumes  on  the  basis  of  optical  property  contrast 
from  measurements  conducted  at  the  air-tissue  interface.  However,  the  ability  to  optically  image  or  detect  diseased 
tissue  volumes  located  deep  within  tissues  depends  upon  the  contrast  provided  by  differences  in  absorption  and 
scattering.  The  exogenous  contrast  offered  by  fluorescent  contrast  agents  may  be  superior  to  that  provided  by 
nonfluorescing,  light-absorbing  compounds,  when  the  optical  measurements  are  conducted  with  frequency-domain 
techniques.  However,  the  reconstruction  of  internal  fluorescent  properties  of  quantum  efficiency  and  lifetime  has  been 
difficult,  especially  when  the  finite  partitioning  of  fluorescent  compounds  takes  place  between  normal  and  diseased 
tissues.  Also,  the  correct  absorption  coefficient  map  is  required  for  the  successful  reconstruction  of  lifetime.  Herein 
we  present  a  novel  fluorescence-enhanced  imaging  algorithm  for  frequency-domain  photon  migration  measurements 
conducted  at  the  air-tissue  interface.  Similar  to  Born  iterative  image  reconstruction  techniques,  fluorescence-enhanced 
imaging  differs  in  that  it  utilizes  measurements  of  generated  fluorescent  wave  instead  of  scattered  excitation  wave. 
Using  synthetic  data  sets,  we  demonstrate  fluorescence-enhanced  imaging  using  FDA  approved  fluorescent  agent. 
Indocynine  Green  (ICG).  Our  results  show  the  fluorescence-enhanced  imaging  algorithm  works  well  up  to  10:1  dye 
uptake  ratio,  and  it  is  relatively  insensitive  to  measurement  noise.  In  addition,  we  present  the  lifetime  reconstruction 
with  a  modified  fluorescence- enhanced  imaging  technique. 

Keywords:  Image  reconstruction,  Fluorescence,  Lifetime 

1.  INTRODUCTION 

Photon  migration,  or  light  scattering  through  turbid  media,  offers  new  possibilities  for  medical  applications  such  as 
biomedical  diagnostics  and  imaging.  In  near-infrared  (NIR)  wavelength  range  (600  nm  -  1300  nm)  the  absorption  of 
light  by  endogenous  chromophores  is  relatively  low,  ranging  from  0.1  ^  1.0  cm”1,  while  the  scattering  coefficient  is 
high  -  typically  100  ~  1000  cm”1.1,2  The  low  absorption  properties  allow  the  transmission  of  measurable  amount  of 
visible  and  near-infrared  light  deep  into  the  tissue  volume.  The  high  tissue  scattering,  however,  reduces  the  contrast 
and  deteriorates  image  resolution  formed  with  this  light,  especially  when  tissue  structures  of  interest  are  buried 
deep  within  the  tissue.  The  endogenous  optical  contrast  for  detection  of  diseased  tissues  can  be  provided  by  (1)  the 
local  absorption  coefficient,  fia ,  which  arises  primarily  due  to  hemoglobin  and  increased  vascular  volumes  associated 
with  tumor  angiogenesis,1  and  (2)  the  isotropic  tissue  scattering  coefficient,  \i\.  The  capabilities  that  improve  the 
detectability  of  the  diseased  tissues  lie  upon  (1)  the  degree  of  optical  contrast,  A  fiQ  and  A/;',  or  the  differences  in 
absorption  and  scattering  properties  that  must  exist  between  normal  and  diseased  tissue  for  the  effective  detection 
and  (2)  the  successful  implementation  ofimage  reconstruction  algorithms. 

In  past  decade,  several  investigators  have  sought  ways  to  reconstruct  internal  tissue  optical  properties  to  differen¬ 
tiate  diseased  tissues  from  normal  tissues  based  upon  the  optical  contrast  due  to  absorption  and  isotropic  scattering 
coefficients.3-5  However.  Troy  ei  oL5  have  shown  from  tbeir  in  vitro  measurements  of  breast  tissue  specimen  that 
the  optical  contrast  due  to  normal  and  diseased  tissues  may  not  be  enough  for  the  successful  detection  of  tumors. 
The  efforts  to  enhance  the  optical  contrast  by  the  use  of  exogenous  contrast  agents  such  as  fluorescent  compounds 
have  received  more  attention  recently.10-12  Sevick  et  a/.13  showed  that  the  exogenous  contrast  offered  by  fluorescent 
compounds  is  superior  to  that  provided  by  nonfluorescing,  light-absorbing  compounds  when  time-dependent  photon 
migration  measurements  are  employed.  While  the  preferential  uptake  of  fluorescent  contrast  agents  into  disease  tissue 
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volume  of  interest  is  responsible  for  contrast,  the  kinetics  of  fluorescence  decay  processes  can  be  environmentally 
specific  to  tissue  volumes  and  further  induce  additional  optical  contrast  for  the  detection. 

In  this  paper,  we  present  a  novel  fluorescence-enhanced  imaging  algorithm  using  frequency-domain  photon  mi- 
eration  measurements.  Using  synthetic  data  sets,  we  demonstrate  the  reconstruction  of  the  absorption  coeffiu  i 
due  to  fluorophores,  and  simultaneous  reconstruction  scheme  for  the  absorption  coefficient  and  llfe^n'e(  ®U™  J 
show  that  the  fluorescence-enhanced  reconstruction  algorithm  works  works  well  up  to  10:1  uptake  ratio,  and 
relatively  insensitive  to  measurement  noise. 


2.  OPTICAL  CONTRAST 

Frequency-domain  photon  migration  imaging  consists  of  launching  sinusoidally  modulated  excitation  light  at  source, 
location,  r.  and  detecting  the  wave  of  excitation  and/or  emission  light  at  detector  position.  ra.  ropagation  o  i- 
incident  excitation  wave  through  the  tissue  results  in  a  decrease  in  amplitude,  M(f,w),  and  an  increase  in  phase 

shift  ,  9{r,u>). 

2.1.  Contrast  due  to  absorbing  heterogeneities 

The  presence  of  a  light-absorbing  heterogeneity  results  in  the  reflection  of  the  propagating  wave  from  its  position .  vh . 
contributing  a  small,  first-order  scattered  wave,  *fcatt(F<i)  to  the  detected  signal.  The  detected  wave  at  the  detector 
position,  u,  is  the  sum  of  the  incident  wave  and  the  scattered  wave,  $fne(»v)  +  $%att(r d) ■  ^ contrast  is  definec 
as  the  detected  signal  in  the  presence  of  the  heterogeneity  minus  that  detected  in  the  absence,  the  contrast  owing  to 
the  absorption  will  be  determined  by  the  scattered  wave,  $featt(Fd)- 

Even  with  the  perfect  absorber  in  the  lossless  scattering  media,  the  dominant  first-order  scattered  wave.  iv). 

is  small  in  comparison  to  the  incident  wave,  Hence  the  contrast  due  to  the  absorbing  heterogeneity  ,s 

expected  to  be  small.  As  the  absorption  contrast  between  the  heterogeneity  and  its  surrounding  decreases,  the 
magnitude  of  the  scattered  wave  becomes  even  smaller.  In  addition,  when  the  heterogeneity  is  buried  deep  uithin 
the  tissue,  the  contrast  available  for  the  successful  image  reconstruction  will  diminish  even  more  and  can  be  lost  ui 
the  noise  of  the  phase  and  amplitude  measurements. 


2.2.  Contrast  due  to  uptake  of  fluorescent  dyes 

Consider  the  perfect  uptake  of  fluorescent  dye  (i.e.  no  exogenous  fluorescence  attributed  to  the  surroundings)  with 
the  assumption  that  the  excitation  and  emission  spectra  are  well  separated.  The  emission  wave  generated  fioni 
the  heterogeneity,  $™et{rd)  can  be  detected  at  the  tissue-air  interface  with  an  appropriate  interference  filter.  The 
emission  wave  generated  by  the  incident  excitation  wave,  $fnc{rh)  is  proportional  to 

where  <?  is  the  quantum  efficiency  of  the  fluorophore;  r,  its  lifetime  and  par-+m,  its  absorption  coefficient  owing  to 
the  fluorophore.  Thus,  we  can  expect  that  amplitude  and  the  phase  lag  of  the  detected  emission  wave  should  change 
according  to  changes  in  <f> ,  and  r.  By  manipulating^,  and  r,  we  could  expect  the  contrast  offered 

by  the  fluorescence  of  the  heterogeneity  to  be  greater  than  that  offered  by  its  absorbance. 

When  the  finite  partitioning  of  a  contrast  agent  occurs  between  diseased  and  normal  tissues,  the  fluorescent  signal 
originated  from  the  surroundings  or  background  has  to  be  accounted  for.  A  series  of  distributed  fluorescent  sources 
also  contribute  to  the  resultant  phase  and  amplitude  measurements  at  the  detector  position  and  affect  the  final 
optical  contrast.  As  the  concentration  of  fluorophore  in  the  surroundings  increases,  the  amplitude  modulation  of  the 
detected  wave  increases,  and  the  phase  lagwilibe  affected  by  the  strong  presence  of  the  fluorophores  in  surroundings. 
However,  in  the  presence  of  a  fluorescent  dye-laden  heterogeneity,  a  second  source  of  fluorescent  generation.  f  (F). 
is  set  up  in  coincidence  with  the  location  of  the  heterogeneity,  r*,  within  the  tissue  phantom.  In  the  presence  of  both 
heterogeneity  and  finitely  partitioned  fluorophores  in  the  surroundings,  the  detected  signal  at  the  detector  consists  of 
$£f(r)  and  As  in  the  contrast  offered  by  absorption,  the  contrast  offered  by  the  fluorescence  is  determined 

bv  the  size  of  the  contribution  of  relative  to 

The  greater  the  partitioning  of  the  fluorescent  dye  within  the  heterogeneity,  the  stronger  will  be  its  fluorescent, 
source  enabling  optimal  contrast.  Fluorescent  dyes  whose  quantum  yield  and  lifetime  changes  upon  the  preferential 
uptake  by  the  diseased  tissue  can  aid  contrast.  For  example,  if  the  dye  experiences  enhanced  and  reduced  partitioning 
measured  by  the  ratio  of  fiax-^m  of  the  heterogeneity  and  its  surroundings,  the  modulation  will  increase  or  decrease, 
respectively.  Even  if  the  quantum  efficiency,  <f>,  and  lifetime,  r,  do  not  change  in  the  heterogeneity  and  in  its 
surroundings,  the  ratio  of  and  the  additional  phase-shift  due  to  lifetime,  r,  will  give  rise  to  the  larger  contrast 

for  measurement  at  the  fluorescent  wavelengths  than  excitation  measurements. 
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3.  FORMULATION  OF  THE  INVERSE  PROBLEM 

In  this  paper,  we  consider  a  simulated  infinite  homogeneous  turbid  medium  with  a  spatially  uniform  distribution  of 
fluorophores  in  the  background.  For  the  purpose  of  this  description,  we  identify  the  normal  tissues  surrounding  the 
diseased  tissue  by  the  ’background’,  while  the  diseased  lesion  will  be  identified  as  the  ’heterogeneity'.  Heterogeneities 
are  simulated  by  regions  of  increased  uptake  and  changed  lifetime.  Excitation  light  is  launched  at  the  air-tissue 
interface  and  propagates  within  the  phantom.  Upon  encountering  a  fiuorophore,  either  in  the  background  or  het¬ 
erogeneity,  the  excitation  photons  activate  the  fiuorophore  and  result  in  the  emission  of  fluorescent  photons  which 
propagate  in  the  media.  We  assume  that  the  dye’s  fluorescence  and  absorption  spectra  are  separated  so  that  we 
can  ignore  the  possibility  of  the  excitation  of  fluorophores  by  the  fluorescent  reemission.  In  addition,  we  ignore  the 
possibility  of  photo  bleaching. 

In  highly  scattering  media  such  as  tissues,  the  propagation  of  light  is  well  described  by  the  diffusion  approximation 
of  the  radiative  transport  equation.14  The  excitation  light  fluence,  <£r(r,w),  or  the  angle  integrated  flux  of  light 
modulated  at  angular  frequency,  u t  at  position  r  is  described  by  the  diffusion  equation  in  the  frequency  domain  for 
excitation  light  (superscript  x): 


V  ■  [D*  (f)V$r(r, w)J  +  H£(f)  +  y)**(r>)  =  -S*(f)  (1) 

where  iixa{r)  is  the  absorption  coefficient,  D*(f)  =  [3  -  (/ix(r)  +  //'/(f)]"1  is  the  optical  diffusion  coefficient,  c  is  the 
speed  of  light  in  the  medium.  Sx(r)  is  the  excitation  source  term,  and  fifsT(r)  is  the  isotropic  scattering  coefficient. 
For  a  point  source,  Sx  =  M*6(r  —  rl),  where  Mx  is  the  complex  source  AC  amplitude,  and  ?\«  is  the  position  of  the 
source. 

The  complex  emission  light  fluence,  <£m(r,Lj),  is  also  described  by  the  diffusion  equation  in  the  frequency  domain. 
The  excitation  fluence  is  coupled  with  the  emission  fluence  as  the  source  term  in  the  emission  diffusion  equation. 
Superscript,  m  denotes  the  emission  characteristics. 


V  ■  +  (-11?  +  y)  <*m(f,w)  = 


where  <p  is  the  quantum  efficiency  of  the  fiuorophore,  fiax^m  is  the  absorption  coefficient  due  to  the  fiuorophore. 
Sm(r, u;)  is  the  emission  fluence  source  term,  and  Dm(r)  is  the  diffusion  coefficient  at  the  emission  wavelength,  and 
r(r)  is  the  probe  lifetime  at  position,  r. 

Since  the  emission  diffusion  equation  is  an  inhomogeneous  differential  equation,  the  Green's  function  is  used 
to  obtain  the  analytical  solution  cf  the  emission  fluence,  $m(r).  The  emission  fluence,  (rd,r9).  at  the  detector 
position,  rd,  and  the  excitation  source  position,  r, ,  is  expressed  as  follows. 


-S^w) 


1  —  iujr(r) 


(2) 


$m(rd,rt)  =  [  GJ(rd,?)Sm{r',rs)dSl,  d£l  =  dr 

Jn 

'(?)( i 


L 


(  Gj(rd,  ?) y  (P,  r~)  dQ 
n  Z)m(r/)(1  —  iu;r(r/)) 


(3) 


where  Q  is  the  integration  space  and  G/ (•)  is  the  Green’s  function  solution  of  the  emission  diffusion  equation. 

In  order  to  reconstruct  the  spatial  map  of  / tax-+m  detailing  the  heterogeneity,  we  discretize  Eqn.(3)  into  a  se¬ 
ries  of  equations  in  terms  of  G^,  and  measurements  of  ^m(r^,r5).  The  solution  of  the  forward  problem  anc! 
determination  of  G/  and  $x,  is  obtained  with  MUDPACK,  a  finite  difference  differential  equation  solver.15.  Both 
forward  and  inverse  problems  are  carried  out  in  a  two-dimensional  geometry,  since  the  three-dimensional  geometry 
^°r  preliminary  computations  presented  herein.  We  simulate  measurements  of  6m  and  Mm  (or 
^  =  M  e  6  )  at  detector  position,  rd,  in  response  to  excitation  source  at  rs. 

By  discretizing  Eqn.(3), 


i=i 


An(fy)(  1  -  iu>r(fj)) 


(4) 
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where  n  is  the  number  of  total  voxels,  h 2  is  the  area  of  one  pixel,  and  i  =  1 . .  .nd,  j  =  1 . .  .Tvoxel,  and  *  -  *  •  ■  n>- 

respectively,  where  is  the  number  of  detectors,  Tvoxel  is  the  total  number  of  voxels,  and  n,  is  the  number  ot 

sources.  The  consequent  matrix  equation,  =  FX,  is 


Fn  •••  Fln  \ 

(  X[fx)  \ 

— 

F2\  •  •  •  F 2n 

X(r2) 

k  $m(r„fi)m  j 

\  Fmi  ■  *  *  Fmn  j 

X(fn)  J 

_  G  f{rdi,rj)$x  {r^J^jr^h2 

ij  ~  ^(fjKl-iwr^)) 

X(rj)  —  Pax— ►m(fj') 

where  F  €  Cmxn,  X  €  Rn ,  €  Cm,  respectively.  Also,  m=n<f  x  nS)  n  is  the  total  number  of  voxels  (Tvoxel=  IT  ) 

In  the  case  of  reconstruction  of  r,  the  matrix,  F2,  and  the  solution  vector,  Y,  are  constructed  such  that 


Yfi) 


Gj(fdi,rj)^£(fsiirj)^ax^m{rj)d{rj)h2 

Dm(fj) 

1 

(1  -  iwr(rj)) 


(6) 


where  F2  €  CmxnJ  Y  €  Cn,  respectively. 


4.  SIMULATION 

The  simulation  of  Forward  and  Inverse  problems  in  lifetime  reconstruction  is  carried  out  using  MUDPACK.  a  finite 
difference  partial  differential  equation  solver.15  Both  forward  and  inverse  problems  are  simulated  on  17x1/  grids, 
and  the  dimension  of  2-D  geometry'  is  chosen  to  be  a  4  x  4  [cm2]  square  leading  to  the  pixel  length  of  0.25  cur. 

The  refractive  index  of  the  phantom  is  set  to  1.33  to  mimic  the  properties  of  tissue.  A  source  is  centered  on  each 
side  of  the  square  phantom,  and  56  detectors  are  evenly  spaced  around  the  phantom.  Under  condition  of  perfect 
uptake,  a  fluorescent  contrast  agent  will  preferentially  be  taken  up  in  a  tissue  volume  of  interest.  However,  we 
confine  our  simulation  to  consider  imperfect  uptake  with  100:1  uptake  ratio  unless  otherwise  specified.  Also,  in 
the  lifetime  reconstruction,  we  will  consider  the  case  where  the  lifetime  of  the  heterogeneity.  rh .  is  shorter  than 
that  of  the  background,  r6.  The  following  table  lists  the  values  of  optical  properties  used  in  the  simulation.  There 


Table  1.  Optical  properties  of  background  and  heterogeneity  in  the  simulation. 


Optical  Properties 

Pam-+t 

(cm-1) 

Pam— ►m 

(cm-1) 

Kr.Km 

(cm-1) 

4> 

(ns) 

Par-41 

(cm'1) 

Par-4  n? 

(cm-1) 

Background 

0.002 

0.0 

10.0 

!  0.03 

1.0 

0.0 

0.002 

Heterogeneity 

0.002 

0.0 

10.0 

j  0.03 

1.0 

0.0  | 

0.2 

are  two  parts  in  the  simulation:  the  forward  and  inverse  simulations.  The  forward  simulator  is  used  to  compute 
simulated  measurements.  For  solutions  of  the  forward  problem,  the  known  optical  properties  of  the  background  and 
the  heterogeneity  shown  in  Table  1  are  used  to  calculate  the  fluence  at  each  node,  and  the  complex  fluence  values 
are  collected  at  the  detector  nodes.  Gaussian  random  noise  is  added  and  the  calculated  values  are  used  as  simulated 
measurement  input  into  the  inverse  problem.  The  following  describes  the  simulation  data  generation  and  inversion. 

♦  Forward  solution 

1.  One  of  four  isotropic  point  sources  of  1/D*(r,)  is  active  at  a  time.  The  isotropic  point  sourqe.  3(pfl.r  4* 
—  r5),  represents  the  incident  light  at  the  simulated  air-tissue  interface. 
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2.  The  complex  excitation  fluence,  <$r(r),  is  calculated  at  each  grid  point  in  response  to  an  isotropic  source. 
Zero  fluence  boundary  condition  is  imposed. 

3.  The  calculated  excitation  fluence,  $r(r),  at  each  node  is  coupled  into  the  source  term  of  the  emission 
diffusion  equation  in  Eqn.(2),  and  the  complex  emission  fluence,  4>m(r, r3),  is  calculated.  Then,  the 
complex  emission  fluence  at  the  detector  position  rrf,  $m(rd,r,)  is  determined  at  56  detector  positions. 

4.  The  next  source  becomes  active  and  the  previous  steps  are  repeated.  In  the  end,  a  total  of  224  simulated 
detector  measurements  are  stored  in  a  column  vector,  $m. 

5.  Random  Gaussian  1%  AC  amplitude  and  5°  (or  10°)  phase  noises  are  added. 

•  Inverse  solution 

The  spatial  map  of  absorption  coefficient  due  to  fluorophores,  /iflr-+m(r)  is  obtained  by  Levenberg-Marquadt 
method.  The  reconstruction  of  lifetime  of  the  fluorophores  follows  the  same  inversion  procedures  except  the 
parameter  of  interest  is  r  instead  of  The  procedure  is  as  follows. 

1.  The  Green’s  function,  G/(,}  and  the  background  excitation  fluence,  are  solved  using  Ml  DPACK.  For 

the  first  iteration,  $x(rs,r)  is  set  equal  to  $f(r,,r)  and  /iax_fm(r)  to  (i  ?-  Born  approxima¬ 

tion).  The  iteration  counter,  k,  is  set  to  zero,  and  the  regularization  parameter,  A.  is  adjusted  to  an  initial 
value. 

2.  The  matrix  is  calculated  from  G/(,),  $x(r,,r),  and  given  estimates  of  <j>,  r,  and  Dm . 

3.  The  update,  A is  obtained  by  solving 

ax  =  (F{k)F(k)H)  +  -  Fik)xik)) 

where  X  is  jjLax-+m,  and,  for  lifetime  measurement,  X  is  equal  to  r.  The  regularization  parameter.  A  is 
incorporated  to  reduce  the  condition  number  of  F^F^  matrix  since  the  matrix  is  ill-posed.10 

4.  (F)  —  JJ-ax-tm  (F)  +  A/iai-4m  (F). 

5.  Calculate  RMSE^l 

6.  If  RMSE^  >  RMSE^'1),  A  is  increased  and  procedures  starting  at  step  3  are  repeated. 

7.  Increment  the  counter,  k=k+l. 

8.  The  forward  solver  is  called  and  the  excitation  fluence,  <$>r(M  is  updated. 

9.  Repeat  until  the  convergence  criteria  is  satisfied. 

This  reconstruction  algorithm  is  different  from  other  conventional  absorption  imaging  modalities  in  that  the  ab¬ 
sorption  coefficient  owing  to  the  fluorophore  is  updated  with  the  detected  fluorescent  wave,  and  not  with 

Thus,  only  the  incident  excitation  wave,  $fnc(r)  is  updated  iteratively,  rendering  the  problem  a  Fredholm 
integral  equation  of  the  first  kind.  This  algorithm  utilizes  the  added  optical  contrast  in  terms  of  phase-shift  and 
amplitude  demodulation  due  to  the  fluorophore. 


5.  RESULTS 

To  compare  the  overall  accuracy  of  the  reconstruction,  the  relative  root- mean-square  error  (RMSE)  of  the  recon¬ 
structed  profile  as  a  function  of  the  iteration  steps.  The  relative  RMSE  is  defined  as 


RMSE{i) 


I'LilXVHfj)  -  X*(rj)]2 


where  X^{fj)  is  the  reconstructed  value  of  the  parameter  of  interest  in  the  ith  iteration  or  r)  and  in  the 

jth  cell,  with  X(rj)  being  the  actual  value.  While  artificial  for  this  synthetic  demonstration  of  the  inverse  imaging 
algorithm,  our  criteria  nonetheless  provides  information  for  comparisons  between  several  test  cases  given  in  following 
sections.  In  the  ultimate  execution.  RMSE  is  determined  from  the  errors  between  measured  and  computed  phase 
and  amplitude  attenuation. 
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Figure  1.  pax->m 
/4r-m  =  0.002  cm- 
comparison  between 


reconstruction  with  1%  amplitude  and  5°  phase  error,  100:1  uptake  ratio.  —  0-2  cm  . 

•i,  frequency=100MHz,  and  A  =  10~8.  (a),  (b)  object  at  (4,4).  and  (c),  (d)  object  at  (9.9).  (e) 
RMSE(i)  values  without  noise  and  with  5°  and  10°  phase  noise  and  1%  amplitude  noise. 


5.1.  Reconstruction  with  noise  _ 

Random  Gaussian  noise  of  1%  amplitude  and  5°  and  10°  phase  errors  are  added  to  4>m  (u)  to  test  the  st  ability  of  the 
algorithm.  Figure  1  shows  that  this  reconstruction  algorithm  was  able  to  reconstruct  the  absorption  map  even  in  the 
presence  of  10°  phase  random  noise.  This  robustness  of  algorithm  is  attributed  to  the  added  phase  and  amplitude 
contrast  from  the  fluorescence  dynamic,  which  is  absent  from  the  conventional  absorption  algorithm.  One  can  see  the 
effect  of  the  added  noise  in  the  RMSE  vs.  iteration  plot.  As  expected,  the  calculated  RMSE  values  in  the  presence 
of  amplitude  and  phase  noises  are  considerably  higher  than  the  case  without  noise.  The  small  difference  in  RM..E 
values  indicate  that  the  presence  of  noise  does  not  affect  the  reconstruction  of  naz-+m  significantly. 

5.2.  r  and  (tax-tm  simultaneous  reconstruction 

In  the  previous  section,  we  showed  that  the  absorption  coefficient  reconstruction  with  the  emission  wave,  4>m. 
instead  of  the  excitation  wave,  was  possible.  Since  the  additional  phase  contrast  comes  from  the  lifetime  of  the 
heterogeneity,  rh,  we  can  consider  the  phase  contrast  originating  from  rh  as  an  additional  phase  error.  Taking  into 
account  that  the  fluorescence-enhanced  imaging  can  withstand  high  phase  error,  we  construct  the  map  of  fiar-*m 
first,  and  then  construct  the  lifetime  map.  The  reconstruction  of  and  r  was  alternately  repeated  until  the 

desired  RMSE  values  were  achieved.  Figure  2  and  Figure  3  show  the  reconstruction  results.  In  the  simultaneous 
reconstruction  scheme,  the  original  Mor-m  value  of  0.2  cm'1  was  recovered  and  r  reconstruction  converged  after  less 
than  20  iterations.  The  final  reconstruction  profiles  indicate  that  it  is  desirable  to  have  an  additional  phase  contrast 
as  well  as  the  absorption  coefficient  change  due  to  the  finite  partitioning  of  fluorophores. 

5.3.  Results  owing  to  perturbation  of  a  strong  absorption  region 

In  real  situations,  there  may  exist  variation  in  the  absorption  coefficient  due  to  the  endogenous  chromophores  inside 
the  phantom.  The  simulation  is  performed  to  test  the  reconstruction  algorithm  in  the  presence  of  a  strong  absorption 
regions  besides  the  heterogeneity  positioned  at  (9,9).  The  reconstruction  results  are  shown  in  Figure  4.  Since  both 
absorption  and  lifetime  reconstruction  is  performed  using  the  emission  wave  information,  the  fluorescence-enhanced 
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Figure  4.  Reconstruction  of  (i4Mm  and  r  in  the  presence  of  strong  absorption  region  at  (4.4).  Heterogeneity  is 
located  at  (9,9).  100:1  uptake  ratio,  the  reconstruction  results  were  achieved  after  200  iterations. 

imaging  algorithm  was  able  to  ieconstruct  both  lifetime  and  absorption  coefficient  maps  due  to  the  preferential  upt  ake 
of  fluorescent  dyes  between  normal  and  diseased  tissue  volumes  despite  the  unaccounted  non-fluorescing  regions  of 
strong  absorption. 


6.  CONCLUSIONS 

In  this  paper,  we  demonstrated  the  novel  image  reconstruction  algorithm  using  additional  phase  and  amplitude  con¬ 
trast  originating  from  fluorescent  contrast  agents.  The  fluorescence-enhanced  absorption  imaging  is  robust  against 
the  phase  and  amplitude  errors,  and  using  this  robustness  of  the  algorithm  against  noise,  we  were  able  to  accurately 
reconstruct  absorption  coefficient  and  lifetime  map  due  to  the  preferential  uptake  of  fluorescent  dyes.  The  simul¬ 
taneous  reconstruction  scheme  converged  to  the  original  values  of  and  r  faster  due  to  the  additional  phase 

and  amplitude  information  in  the  emission  wave.  Moreover,  the  presence  of  the  strong  absorption  regions  did  not 
affect  the  reconstruction  of  lifetime,  which,  in  turn,  indicates  the  robustness  of  this  algorithm.  We  propose  that  the 
development  and  use  of  the  fluorescent  agents  which  change  lifetime  upon  preferential  uptake  into  diseased  tissue 
volumes  will  improve  discrimination  of  diseased  tissue  volumes. 
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Fluorescence  Lifetime  Imaging  and  Spectroscopy 
in  Random  Media 
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Table  4.2.  Optical  Properties  and  Experimental  Parameters  Used  as  Inputs  for  the  Forward 
Problem  Along  with  Values  of  <p  Obtained  from  the  Reconstructions 
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